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'“Several  problems  in  estimation  theory  with  current  application  interest  to 
orbit  determination  are  treated  in  tfciw  dleeertntlon.  Although  much  of  the  material 
:overed  has  general  application,  the  important  specific  of  Interplanetary  orbit  deteml- 
latlon  is  given  special  emphasis  and  is  used  throughout  to  provide  context  and  to  obtalr 
numerical  results.  Topics  covered  Include:  modeling  the  errors  in  the  navigation 
icasurements,  describing  the  navigation  measurements  compactly  in  terms  of  the  unknown 
>arameters,  and  describing  the  estimation  accuracy  in  a  variety  of  poorly  modeled  situ¬ 
ations.  (  sj  — 

The  phase-coherent  doppler  system  used  by  NASA's  Deep  Space  Net  is  modeled  as  a 
time-invariant  linear  system  to  determine  the  influence  the  master  oscillator  'stability 
tas  on  the  resultant  doppler  data  accuracy. 

Analytical  expressions  for  determining  the  errors  Introduced  in  spacecraft  state 
estimates  due  to  Earth-based  tracking  station  location  errors  and  constant  unmodeled 
forces  acting  on  the  probe  are  derived.  This  is  accomplished  through  re-consideration 
ind  extension  of  time  explicit  analytic  expressions  for  the  doppler  and  range  observable 
is  a  function  of  the  probe  state. 

The  Householder  transformation  technique  for  the  solution  to  the  linear  least 
iquares  problem  is  recalled.  This  method  has  been  shown  to  be  ntaerically  superior  to 
Llther  the  classical  normal  equations  or  standard  Kalman  filtering  approaches  to  the 
Lolutlon  of  the  least  squares  problem  by  avoiding  certain  squaring  operations  in  the 
boraputational  algorithm.  Algorithms  are  derived  which  permit  a  wider  range  of  applica- 
ility  of  the  Householder  technique  to  estimation  problems. 
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PREFACE 


The  constant  search  for  the  development  of  efficient  and  improved 
estimation  techniques  for  sensor  data  in  U.S.  Air  Force  Systems  is  an 
effort  which  will  continue  indefinitely  with,  presumably,  ever  improving 
results.  This  report  develops  a  rather  comprehensive  set  of  new  results 
in  this  area  of  fundamental  importance  to  so  many  U.S.  Air  Force  applications 
including  Space  Systems,  Missile  systems,  tactical  systems,  aeronautical 
systems,  and  other  U.S.  Air  Force  systems  although  the  specific  vehicle 
for  application  throughout  this  report  is  space  systems. 

Tills  research  report  was  prepared  under  research  projects  supported 
by  the  U.S.  Air  Force  Office  of  Scientific  Research  under  AFOSR  Grant  72*2166, 
Design  of  Aerospace  Systems,  and  by  the  U.S.  Air  Force  Space  and  Missile 
Systems  Organization  under  Contract  No.  F04701-72-C-0273,  Advanced  Space 
Guidance,  and  this  report  constitutes  part  of  the  final  report  for  these 
contracts . 

The  research  described  in  this  report  "Problems  In  Estimation  Theory 
With  Application  To  Orbit  Determination,"  UCLA-ENG-7275,  by  David  Walter 
Curkendall,  was  carried  out  under  the  direction  of  C.T.  Leondes  and  E.B.  Stear, 
Co-Principal  Investigators  in  the  Schools  of  Engineering  in  the  University  of 
California  at  Los  Angeles  and  Santa  Barbara,  respectively. 
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ABSTRACT 


Several  problems  in  estimation  theory  with  current  application 
interest  to  orbit  determination  arc  treated  in  this  dissertation. 
Although  much  of  the  material  covered  has  general  apolicat  ion,  the 
important  specific  of  interplanetary  orbit  determination  is  given 
special  emphasis  and  is  used  throughout  to  provide  context  and  to 
obtain  numerical  results.  Topics  covered  include:  modeling  the 
errors  in  the  navigation  measurements,  describing  the  navigation 
measurements  compactly  in  terms  of  the  unknown  parameters,  and 
describing  the  estimation  accuracy  in  a  variety  of  poorly  modeled 
situations. 

The  phnse-cohcrcnt  doppler  system  used  by  NASA's  Deep  Space  Net 
is  modeled  os  a  time- invariant  linear  system  to  determine  the  in¬ 
fluence  the  master  oscillator  stability  has  on  the  resultant  doppler 

d'.tfs  r.ecurvay. 

Analytical  expressions  for  determining  the  errors  introduced  in 
spacecraft  state  estimates  due  to  Earth-based  tracking  station  loca¬ 
tion  errors  and  constant  unmodeled  forces  acting  on  the  probe  an 
derived.  This  is  accomplished  through  re-consideration  and  extension 
of  tine  explicit  analytic  expressions  for  the  doppler  and  range  ob¬ 
servables  as  a  function  of  the  probe  state. 

lhs  Householder  transformation  technique  for  the  solution  to  the 
linear  least  squares  problem  is  recalled.  ThJi  method  has  bean  shown 
to  be  numerically  superior  to  cither  the  classical  normal  equations  or 
standard  Kalmar  filtering  approerhrs  to  the  solution  of  the  least 
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squares  problem  by  avoiding  certain  squaring  operations  in  the  compu¬ 
tational  algorithm.  Algorithms  are  derived  which  permit  a  wider  range 
of  applicability  of  the  Householder  technique  to  estimation  problems. 
Specifically,  methods  are  provided  for  calculating  the  actual  estima¬ 
tion  performance  when  the  apriori  modeling  used  by  the  filter  is 
incorrect  or  certain  stochastic  processes  affecting  either  the  data 
directly  or  the  probe  state  itself  are  present.  These  algorithms 
are  used  in  several  key' numerical  examples  calculating  the  effects 
of  stochastic  station  location  motion  and  probe  state  i»rocess  noise. 

The  problem  of  detecting  the  presence  of  parameters  not  modeled 
by  the  filter  through  the  inspection  of  the  data  residuals  is  consid¬ 
ered.  Again,  the  Householder  transformation  technique  is  used,  this 
time  to  make  the  problem  of  computing  the  expected  sum-of-squarcs  of 
the  residuals  ns  a  function  of  t!  e  unmodclcd  parameters  computation¬ 
ally  tractable.  Several  suggestions  arc  made  concerning  this  detec¬ 
tion  problem,  and  an  !nfor~?l  i:yv>thcr!s  ierti'p  rr'ccdu”?*  f~r 
validating  least-squares  solutions  is  detailed. 
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CHAPTER  I 


INTRODUCTION 

Navigation  is  an  integral  part  of  the  design  and  execution  of 
any  space  mission.  The  ability  to  determine  and  control  the  flight 
path  of  the  spacecraft  to  high  accuracy  is  a  prime  mission  require¬ 
ment  in  order  that  the  mission  objectives  can  be  met.  In  all  but 
the  simplest  unmanned  missions,  orbit  determination  and  trajectory 
control  are  needed,  not  just  during  the  initial  powered  flight  phase, 
but  throughout  the  mission,  so  that  continuing  refinements  in  know¬ 
ledge  or  changes  to  the  flight  path  can  be  made  to  optimize  the 

scientific  return.  For  manned  missions,  the  additional  concerns  for 

•  r  ;  j 

crew  safety  add  to  the  imperative  reasons  for  these  functions. 

The  current  state  of  the  art  of  navigation  systems  design  re- 
suit  in  mission  performance  limited  not  by  the  mechanization  errors 
in  the  sub-sys terns  that  control  the  changes  in  the  flight  path,  but 
by  the  orbit  determination  accuracy.  The  above  statement  is  a  gen¬ 
eralization  and  cannot  be  taken  categorically  since  it  is  possible  to 

i  1 1  .  i  e  .  . 

i  i'  .  ; 

construct  sensible  mission  situations  where  the  reverse  would  be 

.  •  .i 

true.  Still  in  all  demanding  Interplanetary  missions  currently 

.  * 

/  “  <  J  A  '  *  *  1  .  \ 

planned  -  Mariner/tors  ’71|  Kariner/Venus/Mercury  '73»  Viking  ’75» 
Grand  Tour  Missions,  etc.  -  a  factor  of  10  improvement  in  the  control 

4  ■  '  '  ■  ■■  '  V  *  v  VO  .  -  .j'.  ,  ;;  .  .  v 

sub-systems  would  refine  mission  performance  while  a  similar  factor 

*  t  *  ■■  ^  k  * '  *.  ii  ♦  .  *  » 

*  *-*  f  J)  -•*  *V/ 

in  orbit  determination  accuracy  would  revolutionize  it.  A  more  com- 
plete  description  of  the  state  of  the  art  nnd  current  inherent 
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limitations  of  planetary  navigation  la  found  in  Reference  1. 

Intimation  theory  as  applied  to  orbit  determination  la  selected 
as  the  subject  for  the  dissertation  beeanae  of  this  Importance  to 
overall  mission  design  and  performance.  Several  key  problems  rela¬ 
tive  to  earth-based,  radio  navigation  are  Investigated  and  nee  re¬ 
sults  are  obtained  from  both  the  estimation  theoretic  and  physical 
model  points  of  vise.  A  discussion  of  the  content  and  major  results 
contained  In  the  following  chapters  Is: 

Chapter  II 

the  phase-coherent  doppler  tracking  system  la  modeled  as  a 
linear,  time-invariant  system  and  the  effects  that  variations  in 
the  output  frequency  of  the  master  oscillator  have  on  the  doppler 
measurement  are  determined.  It  la  found  that  the  slow  drifts  in 
the  oscillator  frequency  are  of  greater  Importance  than  the  rela¬ 
tively  much  larger,  but  more  rapid,  frequency  excursions  arising 
from  the  vide -band  component  of  the  oscillator's  output  spectrum. 

In  general,  the  results  presented  perait  the  determination  of 
the  oscillator  performance  characteristics  required  to  achieve  a 
specified  doppler  data  accuracy. 

Chapter  III 

8erles  expansion  methods  are  used  In  order  to  obtain  a  time- 
explicit  analytic  relationship  for  the  doppler  observable  In  terns 
of  the  (unknown)  probe-state  and  the  earth-based  tracking  station 
coordinates.  This  la  sn  extension  of  earlier  work  by  Hamilton  and 
Melbourne  (Reference  2).  The  resulting  expansion  is  used  as  the 
baals  for  on  analytic  orbit  determination  accuracy  analysis  which 
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reveals  the  structure  of  the  accuracies  to  be  expected  as  a  function 
of  the  probe -Earth-8un  geometry.  The  analysis  is  approximate  and 
only  valid  for  relatively  short  spans  of  tracking  data  (less  than 
one  week  for  a  probe  in  interplanetary  space)  but  is  highly  useful 
in  gaining  intuitive  understanding,  establishing  the  key  error 
sources,  and  interpreting  resulta  from  more  precise  analyses.  The 
effects  of  augmenting  the  basic  doppler  data  type  vlth  ranging  data 
are  explored. 

Chapter  IV 

The  analysis  of  the  previous  chapter  is  employed  to  derive 
analytic  expressions  for  the  errors  in  the  estimate  of  probe  state 
as  a  function  of  the  two  most  Important  model  errors:  l)  station 
location  errors  and  2)  unmodeled  forces  affecting  the  probe's 
motion.  This  analysis  assumes  doppler-only  tracking;  similar  ana¬ 
lytic  results  for  the  case  where  ranging  data  are  also  employed  do 
not  appear  tractable.  Because  the  structure  of  the  analytic  analysis 
suggests  that  ranging  data  could  reduce  the  error  sensitivities  to 
the  unmodeled  forces,  a  numerical  study  is  presented  to  test  this 
hypothesis.  The  results  shorn  that  the  position  errors  are  reduced 
by  this  additional  data  in  the  vicinity  of  the  dots  itself,  but  at 
the  expense  of  larger  velocity  errors.  These  results  lay  the 
foundation  for  the  intelligent  application  of  this  supplementary 
data  type  in  operational  situations. 

Chapter  V 

A  short  history  of  the  development  of  square-root  filtering 
using  Householder  transformations  is  r. iv?n  r. r.d  the  applicability  of 
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this  sore  numerically  a  tab  It  computational  approach  la  extended  tc 
lncludt  two  important  problem.  They  art:  to  evaluate  true  fllttr 
performance  vhtn  l)  tht  apriorl  uncertainty  asalgntd  to  the  esti¬ 
mated  variable  in  filter  initialisation  la  Modeled  improperly  and 
2)  Markovian,  stochastic  processes,  of  which  the  filter  is  unaware, 
are  effecting  the  data  and/or  the  state.  These  extensions,  when 
combined  with  algorithm  currently  available,  provide  the  basis  for 
the  investigation  of  the  succeeding  chapter. 

Chapter  VI 

The  limitations  of  the  analytic  analyses  of  Chapters  XXX  and  IV 
permit  computations  involving  neither  long  spam  of  tracking  data 
nor  the  evaluation  of  the  effects  of  stochastic  processes,  to  cir¬ 
cumvent  this,  a  more  precise,  but  non-analytlc,  model  is  adopted  to 
Investigate  the  effects  of  stochastically  varying  station  location 
errors  and  unmodeled  process  noise  as  a  function  of  the  length  of 
the  data  arc.  It  is  found  that: 

1.  Increasing  the  length  of  the  data  arc  decremes 
the  errors  due  to  both  fixed  and  random  station 
location  errors  in  two  of  the  three  components  of 
position.  The  error  in  the  third  coordinate  (probe 
declination)  is  remarkably  insensitive  to  data  span 
and  shows' little  variation  over  the  entire  range 
Investigated  (up  to  k  months  tracking). 

2.  Solving  for  the  fixed  components  of  station  location 
errors  produces  better  performance  even  in  the  pres¬ 
ence  of  the  rrndom  components  which  are  Ignored  by 
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the  filter.  But  again,  the  improvement  is  selective 
(e.g.,  a  factor  of  2.3  in  probe  right  ascension  but 
only  a  factor  of  1.3  in  declination). 


3.  The  sensitivity  of  estimation  error  to  constant, 
unmodeled  forces  increases  by  a  factor  of  four  from 

•f 

that  predicted  by  the  short-term  analytic  analysis 
after  approximately  three  weeks  of  data  are  processed. 
However,  qualitative  structure  of  the  analytic  result 
in  terms  of  which  components  of  the  solution  are 
affected  is  found  to  be  accurate  even  for  relatively 
long  arcs  of  data. 

4.  For  a  given  magnitude  of  the  acceleration  affecting 
the  prebe,  unmodeled  stochastically  varying  accelera¬ 
tions  produce  approximately  a  factor  of  from  ten  to 
forty  greater  estimation  errors  than  do  unmodeled 
constant  accelerations.  However,  by  including 
ranging  data  and  extending  the  filter  to  Include 

.solutions  for  constant  accelerations,  the  effects 
of  the  stochastic  forces  can  be  reduced  dramatically 
(by  a  factor  of  ten  in  the  case  studied.) 

Because  the  above  results  are  not  obtained  with  an  analytic 
analysis,  it  is  difficult  to  determine  their  universality  in  diff¬ 
erent  mission  applications  with  different  geometries.  However,  it 
is  thought  that  the  techniques  developed  are  applicable  to  a  wide 
range  of  interplanetary  orbit  detersd  nation  problems  and  can  be  a 
significant  aid  to  forging  r.pjyoprinte  estimation  strategies. 
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Chapter  VII 

Householder  transformation  techniques  are  used  to  determine  the 
behavior  of  the  minimum  weighted  sum  of  squares  of  ths  data  residuals 
in  ths  presence  of  unoodeled  parameters.  In  sddition,  formulae  are 
derived  which  yield  ths  decrease  in  ths  sum-of-squares  which  would 
be  obtained  by  including  additional  parameters  in  the  solution  with¬ 
out  actually  performing  the  extended  solutions.  Thsse  results  are 
employed  in  a  suggested  informal  hypothesis  testing  procedure  whose 
aims  are  to  arrive  at  a  reliable  data  det  in  ths  sense  that  inspec¬ 
tion  of  the  residual  behavior  trill  reveal  modeling  difficulties  and 
to  conduct  the  search  for  the  source  of  the  difficulty  efficiently. 
Appendix 

Much  of  the  dissertation  assumes  a  general  familiarity  with  the 
earth-based  data  types  being  analysed  and  estimation  methods  used  in 
orbit  determination.  This  background  material  is  reviewed  in  the 
Appendix.  Both  range  and  doppler  data  are  discussed  and  a  brief 
presentation  of  how  the  orbit  determination  problem  is  structured 
for  computational  solution  using  digital  computers  is  given.  Simple 
derivations  of  classical  weighted  least  squares  snd  sequential 
filtering  are  presented  and  these  methods  are  related  to  the  square- 
root  filtering  algorithms  employed  throughout  the  dissertation. 
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CHAPTER  n 


THE  INFLUENCE  OF  MASTER  OSCILLATOR  STABILITY  OR  DOPPLER  DATA  ACCURACY 

« 

The  main  data  type  obtained  during  Earth-baaed  spacecraft  track¬ 
ing  ia  the  two-way  coherent  doppler  aignal  which  ia  proportional  to 
the  relative  range  rate  between  the  atation  and  the  a  pace  craft.  In 
actual  practice,  this  doppler  tone  ia  integrated  by  aeana  of  an  elec¬ 
tronic  counter,  producing  the  raw  data  type  which  ia  proportional  to 
the  rar^e  change  between  these  two  parti  cipanta  over  the  counting  in¬ 
terval.  A  aignificant  contributor  to  the  total  doppler  ayatea  error 
ia  the  frequency  instability  in  the  atation  neater  oscillator  which 
aupplica  the  frequency  reference  for  both  the  original  tranaalssion  to 
the  spacecraft  and  for  comparison  with  the  return  signal  to  detect  the 
doppler  tone. 

The  objective  of  this  chapter  is  to  apply  linear  ayatea  rnalysis 
techniques  to  the  doppler  ayatea  in  order  to  dete  mine  tho  quantita¬ 
tive  relationship  between  the  statistical  description  of  the  oscilla¬ 
tor  performance  and  the  resultant  integrated  doppler  data  errors. 

Both  two-way  coherent  tracking  (Just  described)  and  three-way  tracking 
(where  the  receiving  station  la  at  a  distant  location  and  by  necessity 
eaploys  an  independent  reference  oscillator  for  the  doppler  detection) 
are  investigated.  The  case  of  white  noise  in  the  oscillator  output  is 
treated  quite  thoroughly,  but  it  ia  shown  that  the  resulting  data 
noise,  while  quite  visible  in  the  data  is  not  as  serious  as  that 
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arising  from  an  oscillator  error  output  of  primarily  low  frequency 
content.* 

Linear  System  Analysis 

From  the  dopplcr  system  diagram  shown  and  discussed  in  the 
Appendix,  the  two-way,  coherent  tracking  system  can  be  modeled  as 
shown  in  Figure  2.1.  Here,  the  error  signal,  f^(t)  ir  fed  into  a 
delay  line  of  duration  4  (the  round  trip  light  time  to  the  probe)  and 
simultaneously  fed  toward  to  be  differenced  with  the  output  of  the 
delay  line.  The  result  is  the  instantaneous  doppler  error,  p(,  which 
in  turn  is  integrated  to  produce  a  range-change,  or  integrated  doppler 
error,  AP(.  The  lower  diagram  (ii)  is  the  same  system  and  constitutes 
the  definitions  of  h^(t)  and  hg(t),  the  unit  imnulse  response  of  the 
two  subsystems.  From  inspection: 

h^t)  -  fi(t)  -  6(t  -  i)  (2.1) 

Assume  the  application  of  a  mean  zero  vide  sensc-(or  covariance-) 
stationary  process  of  spectral  composition,  Gf(«)(  the  mean  zero 
assumption  is  not  really  necessary  since  it  is  readily  seen  that  the 
steady-state  output  of  the  system  to  a  constant  input  is  zero).  Re¬ 
calling  the  relation 

0.(«)  -  Gf(«)  |1IA  (j*)|2 

*i#e  "have  to  contend  here  with  the  rather  confusing  semantic  problem 
that  the  error  output  is  a  frequency.  In  modeling  this  error  as  a 
stochastic  process  vc  frequently  must  refer  to  its  power  spectral  den¬ 
sity  and  speak  of  its  frequency  content.  This  dual  use  of  the  term 
frequency  is  unfortunate  but  apparently  unavoidable.  We  shall  attempt 
however,  to  speak  particularly  carefully  in  this  regard  to  keep  the 
confusion  at  a  minimum. 
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where  H^(ju)  ia  the  aysten  function  and  ia  defined  aa  the  Fourier 
transform  of  h^(t). 

-  (  [®(t)  -  6(W)]e’Jtt*dt 

. .  -a- 

IH^dw)!2  -  U«) 

■  2  (1-cot  out 

Therefore,  if  in  particular,  Gf(u))  ia  represented  by  a  white  noise 

source  of  constant  spectral  content,  K 

Q.(cu)  »  2  K  (1-cos  o>4)  (2.2) 

P 

and  would  appear  aa  shown  in  Figure  2.2.  Of  particular  note  ia  the 
attentuation  made  of  the  low  frequency  content  of  the  input  signal. 
This  ia  particularly  important  because  the  second  half  of  the  system 
is  a  pure  integrator,  an  unstable  element.  Thus, 


h. 

1 


*-»> 

Figure  2  2.  Instantaneous  Donpier  Signal  Error  From  a  White  Noise 
Oscillator  Error  Model. 
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suppose  that  G^(cu)  contained  a  component  of  (real  oscillators 

are  supposed  to  contain  a  component  —  ,  the  so-called  "flicker- 

« 

noise"  component).  Since 

2  2 

l-COS  oii  m  ®£“- 

the  spectral  height  of  G.(o)  would  still  be  finite.  This  observation 

P 

is  in  the  nature  of  an  aside  and  may  appear  unmotivated,  but  it  helps 
to  demonstrate  the  important  point  that  bias  errors  in  the  doppler 
dsta  are  extremely  low  in  magnitude.  To  demonstrate  their  existence 
at  all,  second  order  arguments  must  be  pursued. 

Continuing  with  the  white  noise  postulate,  it  is  useful  to  take 
the  inverse  Fourier  transform  of  (2)  to  obtain  the  autocorrelation 
function  of  the  doppler,  i.e. 

R.(t)  -  «u 


"  f U-e'^-e^e^d  » 

■![«”•  I 


■  2  *  (T  ■  *)  ’  ^  6  (t  ♦  j e) 


J  (2.3) 


This  function  is  normalized  and  displayed  in  Figure  2.3. 
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T’ - - - 


c 


The  power  spectrum  of  the  output  of  the  entire  system,  i.«.,  the 
spectrum  of  the  range-change  or  integrated  doppler  errors,  is  com¬ 
puted  as  follows 

<*•>  ‘  3T 

^(Ju.)  -  (1  -  e  m*°l)/fa 

I^O)!2  -  2  (1-cos  )/u>2 

Thus 

G  (u>)  -  2K  (1-cos  «uX)/w2  (2.4) 

dp 

Again,  it  is  useful  to  inspect  the  behavior  near  u»0.  Employing 

L'Hospital's  rule 

d( 1-cos  m£)  ^  As in  ml 

d(«,2)  2«> 

Re-employing  on  the  first  result, 

dUsin  u)Xj  t2 cos 

d(2uj)  2 

this  last  result  becomes  i2/2  at  ^-0.  Thus 

Ga  (u>)  -  K*2 

Ap 

<u-0 

Equation  (4)  is  plotted  in  Figure  2.4.  The  plot  is  normalized 
by  setting  K*l/2  and  A*l.  Then  to  show  the  effect  of  increasing 
station-probe  round-trip  light  time,  the  curve  for  2  is  also 
plotted.  The  essential  feature  here  is  that  the  low  frequency  content 
is  sharply  heightened  with  increasing  l. 

Thus  far,  only  processes  in  the  steady  state  have  been  considered; 
the  theorem  which  states  that  covariance  stationary  random  processes 
operated  on  by  stable,  linear,  time  invariant  systems  will  also  pro¬ 
duce  (in  the  steady-state)  covariance  stationary  processes  has 
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implicitly  been  used.  Recognising  that  tracking  systems  operate  in 
disjoint  time  segments  of  finite  length,  it  is  necessary  to  derive 
statistics  on  ^pe  including  transient  behavior.  Assume,  then,  that 
the  doppler  portion  of  the  system  has  been  operating  for  some  time 
and  has  achieved  steady-state.  At  t*o,  a  switch  is  thrown,  activating 
the  integrator,  and  measurements  commence  (see  Figure  2.5). 


R,  (V  t2) 


yv  v 


Figure  2.5.  Stationary  Doppler  Process  is  Applied  to  Integrator 
at  t«0. 

Following  an  approach  suggested  by  Papoulis  (Reference  3,  Page 
311)  for  obtaining  the  autocorrelation  function  of  a  process  passed 
through  a  linear  system,  first  consider  an  abstract  system  as  shown 


in  Figure  2.6. 


x(t) 


y(t)  -  Lt(x(t)) 

Figure  2.6.  Linenr  System  with  Input  x(t). 

Computing  the  c*oss-covarlance, 
x(tx)  y(t)  -  Lt(x(tx)  x(t)) 

E^tj)  y(t)J  £  nxy(tl’  t2)  "  E  Mx(tl)x(t)J 

interchanging  the  order  of  the  -xpcctaMon  and  the  system  operators, 


the  following  intermediate  result  is  obtained: 

yv)  -  vw 


(2.5) 


Ih 


Leaving  the  abstract  notation  and  returning  to  the  physical  prob¬ 
lem,  note  that  the  L.  is  the  integral  operator  and  R  ia  the  lnstan- 
taneoua  doppler  error  autocorrelation  function  given  by  (2.3).  Thus, 

RpAp(ti»  V  "  l  f  2[«  <Vt) 

o 

8;4p(ti’  *2>  ■;  {0(VH>  u<  vx)  *V<V'» 

-i  (2.6) 

0  t<0 

where  U(t)  ■ 

1  t*0 

the  unit  step  function. 

Equation  (2.6)  is  an  intermediate  result;  the  autocorrelation  of  the 
Integrated  output,  R&p  (t^,  tg),  is  desired. 

Reverting  to  the  abstract  system  notation  and  continuing  from 
Equation  (2.3): 

VW  •  ®  jr{ta)  ] 

-  E  j"  Ltx(t)  y(t2)J 

Again  interchanging  the  operators 

*^y(t,  t2)  -  Vxy(t»  t2)  (2'7) 

Or,  in  the  notation  of  the  physical  system  and  substituting  the  in¬ 
termediate  result  (2.6)  into  (2.7) » 

*2*  ■  if  1  »<  VtJ  -i  «(<•-'.)  U(t2-(t-l))-i  v(t2-(t+j))dt 

°  (2.8) 

r  f1 

The  evaluation  of  (2.8)  is  rather  tedious,  but  straightforward. 

Omiting  the  algebraic  detail,  the  final  result  becomes: 

■*  *  **  •  1  ,  '■s*  ,  .  *  , 

R&0(ti»  V  { *2  *  u*Vhl 
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This  is  the  desired  result:  the  autocorrelation  of  the  integrated 
doppler  error  assuming  a  stationary  input  to  the  final  Integrator  and 
■ensuring  from  t-0.  Equation  (2.9)  is  difficult  to  regard  directly 
and  arrive  at  any  sensible  interpretation.  If,  however,  we  set 


a  simple  result  for  the  variance  of  the  process  is  obtained 

R  .  (t,  t)  -  -/t-(w)U(t-*)\ 

ApAp  tr  i  ) 


or 

*  ,  <t,  t)  -  - 1  t  <  i  (2.10) 

ApAp  tr 


Equation (2. 10)  is  plotted  in  Figure  2.7. 
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the  autocorrelation  rises  until  tj-i  and  then  remains  constant  at  .5 
until  t^  nears  tg.  It  then  rises  sharply  to  obtain  a  maximum  value  of 
1  (which  checks  with  (2.10))and  then  falls  to  a  constant  value  of  .5 
for  all  t^  a  tg  ♦  i.  '  We  are  now  in  a  position  to  sake  and  demonstrate 
two  assertions: 

1.  Por  t^  a  *,  apc  (t)  is  a  wide  sense- (or  covariance-)  sta¬ 
tionary  process. 

In  order  for  the  assertion  to  be  true,  the  process  Ap((t)  for  £ 

«  to*  eust  have  the  following  properties: 


i. 

|E[4p,(t)J 

11. 

K*,<h)  w}* 

ill. 

^,(0]  - 

a  constant 

iv. 

-  ■ 

)(i,,(t2)  .  .)] 

■  Y  (*!»  tg)  - 

V  (|»!  -  tg|) 

Properties  (i)  -  (iii)  clearly  obtain  from  what  has  already  been  said 
and  from  the  fact  that  since  the  input  process  has  mean  zero  and  the 
initial  conditions  are  zero,  the  output  eust  have  mean  zero  also. 
Property  (iv)  can  be  checked  by  re-writing  the  general  autocovariance 
expression  (2.9)  using  the  constraint  t^,  tg  >  4.  The  following  re¬ 
sult  obtains: 

H*  nF  ou- w}  <a-u> 

where  tj,  tg  *  £ 

and  t  -  tA  -  tg 

Thus  £pc(t)  is  a  wide-sense  stationary  process. 
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Thus  the  variance  goes  linearity  vith  tine  until  t  equals  the 
round-trip  light  tine  (or  delay,  l),  at  which  point  the  steady-state 
la  achieved.  This  systen  has  the  rather  unusual  feature  of  achieving 
steady-state  in  finite  time.  With  this  Insight,  return  to  the  general 
expression  (2.9),  hold  t g  fixed  and  inspect  the  behavior  versus  t^. 

A  particular  example  is  shown  in  Figure  2.8.  Here  K  ■  n,  1  ■  1  for 
normalization.  As  shown 


Figure  2.8.  RA0  (t^,  tg)  Versus  t^ 


2.  for  t^,  tg  2  it  Ap(  can  be  m tuned  to  be  generated  by  two 
Independent  processes;  i.e.  Ap((t)  can  be  written  as: 

AP<(t)  -  a  ♦  b(t)  (2.12) 

where  a  is  a  Man- zero  random  variable  of  variance 

4  •  51  (2-u) 

but  acts  as  a  bias  on  any  one  experiment,  and  b(t)  is  a  mean  zero  wide 


sense  stationary  stochastic  process,  independent  of  a  and  having  auto* 


correlation  function. 

iigJlO  u  (A.|T[)  (2.14) 

The  validity  of  this  assertion  is  shown  by  using  (2*13)  and  (2.14)  in 
(2.12)  and  demonstrating  the  original  correlation  function  (2.11)  is 
obtained.  Equation  \2.lU)  is  shown  plotted  in  Figure  (2.9) 


”1  T  -«  A 


Figure  2.9*  Autocorrelation  of  One  Component  of  R  (t) 

dp 

Further,  a  stochastic  process  having  Identical  first  and  second 
moment  characteristics  to  b(t)  can  be  generated  simply  by  passing 


unite  noise  through  a  smoother; 


i.e. 


v(a)dor 


(2.15) 


where  w(u)  is  a  white  noine  process  of  unity  power  spectral  density. 

In  summary, 

1.  the  output  of  the  integrated  dopplcr  system  due  to  white 
master  oscillator  frequency  variations  is  described  in  general  by  a 
non- stationary  process,  but  this  process  reaches  a  steady-state  in 
finite  time. 

2.  If  the  output  is  sampled  only  after  steady-state  has  been 
achieved  (|  seconds  after  initiation),  the  error  can  be  described  for 
visualization  or  perhaps  simulation  purposes  as  having  been  generated 
by  a  bias  (random  between  experiments  and  with  mean  zero,  but  constant 
during  a  single  experiment)  and  a  white  noise  smoother  as  described  by 
(2.15).  Each  of  these  sources  contribute  equally  to  the  variance. 

The  variance  of  the  process  is  directly  proportional  to  the  round- 
trip  light  time. 

Three -way  Data 

Thus  far,  only  the  case  vherc  a  single  trucking  station  both 
broadcasts  and  receives  the  signal  has  been  considered.  Often,  a 
second  station  receives  the  signal  os  it  is  broadcast  from  the  space¬ 
craft.  This  is  done  usually  to  provide  better  geometry  for  near- 

4 

Earth  tracking  (less  then  one  lunar  radius)  since  the  stations  can  be 
separated  by  as  much  as  an  Earth  diameter.  This  configuration  is  the 
so-called  three-way  data  taking  mode,  because  there  are  three  dis¬ 
tinct  participants.  This  new  system  can  be  modeled  in  a  manner  analo¬ 
gous  to  Figure  2.1.  This  is  shown  as  Figure  2.10.  Here,  independent 
oscillators  are  nc-ded. 
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Receive 

Figure  2.10.  Block  Diagram  of  Three-Way  Doppler  System 

and  there  can  be  no  feed  forward  of  the  input  errors  to  ultimately  can¬ 
cel  and  keep  the  output  a  stationary  process;  the  aystem  is  now  unstable. 
In  addition,  the  original  assumption  that  f((t)  contained  no  bias  or 
offset  errors  cannot  be  so  blithely  made;  for  three-way  data,  an  os¬ 
cillator  bias  directly  produces  a  bias  in  the  doppler,  p. 

Even  so,  suppose  the  estimation  process  tc  be  not  sensitive  to 
the  bias  and  only  the  stochastic  errors  in  the  data  to  be  of  real 
concern.  It  is  desired  then,  to  calculate  the  effects  of  a  white 
noise  oscillator  input  on  £p  as  was  done  before.  The  calculations 
proceed  similarly  but  are  much  simpler;  therefore  only  the  results  are 
stated.  If  both  oscillators  have  constant  spectral  density,  K,  the 
autocorrelation  function  on  the  doppler  is  simply 

R.(*0  -  *  ft(r)  (2.16) 

P  i» 

Thus,  if  the  output  were  to  be  smoothed  over  short  time  intervals 
(<4)  so  that  we  may  meaningfully  speak  of  variance,  the  variance  on 
the  three-way  data  is  identical  to  that  of  its  two-way  counterpart 
(compare  with(2.3)). 

Often,  analyser  of  doppler  errors  proceed  no  further  than  this 
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since  for  conceptual  purposes  the  data  can  be  considered  as  direct  p 
measurements  by  suitably  sampling  and  differencing  the  output  of  the 
integrator.  The  importance  of  the  negative  correlation  spikes  in 
(2.3)  can  only  be  underlined,  however,  by  passing  the  three-way  data 
through  the  final  integration  stage  and  comparing  the  results  there. 

Clearly  (2.16)  indicates  that  the  pe  is  a  simple  white  noise 
process,  which  when  integrated  yields  a  Wiener  process,  i.e: 


R 

dp 


*1  <  *2 


*2**1 


or  the  variance 

R  (t,t)  -  -  t  for  all  t^O 

dp  tt 

Thus  while  the  two-way  error  is  limited  by  the  round-trip  light  time,  1, 
the  three-way  error  variance  continues  to  build  throughout  the  track¬ 
ing  interval.  For  a  typical  light  time  of  a  1000  sec.  and  a  tracking 
Interval  of  twelve  hours,  the  relative  variance  of  the  oscillator  noise 
on  three-way  as  to  two-way  data  would  be  in  a  ratio  of  approximately 
30  to  1. 


Numerical  Results 

Doppler  tracking  data  was  obtained  from  the  Pioneer  VI  spacecraft 
at  one  second  sampling  rates.  Successive  samples  were  differenced  and 
divided  by  the  count  time  (in  this  case  unity)  to  produce  smoothed  ‘p 
samples.  The  autocorrelation  function  of  these  data  was  taken  viith 
results  as  shown  in  Figure  2.11. 


0  10  20  30  40  90  40  70  10  90  100  110  120 
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Figure  2.11.  Normalized  Autocorrelation  Function  for  Pioneer  VI 
Doppler  Tracking  Data. 

A  negative  spike  is  expected  in  the  function  at  a  lag  near  the 
round-trip  light  time,  which  in  this  case  was  21.3  *ec.  as  predicted 
by  the  theory  (c.f.  Equation  (2.3))*  This  expectation  is  clearly  re¬ 
warded  as  seen  in  Figure  2.11.  The  negative  spike  at  one  lag  time  is 
due  to  another  phenomenon,  receiver  phase  Jitter,  and  can  also  be  ex¬ 
plained  by  a  simple  model  (see  Curkendall,  Reference  4),  but  is  of  no 
consequence  to  the  present  pursuit.  It  does  serve  to  separste  the  re¬ 
ceiver  Jitter  contribution  to  the  high  frequency  noise  from  that  aris¬ 
ing  from  the  oscillator.  This  aeparotion  yields  that 

-  -  (.03U)2  hz2*  (-  is  defined  by  (2.3)) 

It  It 

♦Here  a  cycle  is  one  cycle  of  integrated  dopplcr  which  at  the  S  band 
broadcast  frequency  of  2200  x  10**  hz  is  equivalent  to  6.5  cm.  in  range 
change. 
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This  would  arise  from  an  rms  deviation  in  the  master  oscillator, 

of  1.4  x  lO-11  for  a  1  sec.  averaging  time.  Other  similar  experi- 
nt 

ments  yielded  deviations  as  high  as 

-  -  (.045)2  h*2 

n 

or  2x  lO"**  jjj-  1  sec.  smoothing  time. 

Employing  this  later  figure  in  the  integrated  doppler  variance 
(Equation  2.10), 

RA0  -  (.045)21 

assuming  1  to  be  20  min.  and  converting  to  metric  units, 

oA  ^  VRApUtt)  -  (.045)  7*=  X  V1200  »ec  X  6.5  ~ 

ysec  y 

■  10  cm. 

The  total  Integrated  doppler  error  from  all  sources  (see  the  discus¬ 
sion  in  ..ppendix  II)  is  approximately  150  cm;  thus  the  white  noise 
component  from  rubidium  oscillators  contributes  but  a  small  part  of 
this  total  error.  But  the  longer  term  drifts  also  characteristic  of 

these  same  oscillators  can  be  shown  to  be  a  major  contributor.  A 

12 

reasonable  drift  over  a  12  hr  pass  is  5  parts  in  10  ,  much  smaller 
than  the  measured  short-term  fluctuations  of  2  parts  in  1011.  Its 
effect  can  be  computed  in  the  following  manner: 

Observe  that  Ap(  can  be  written  directly  in  terms  of  the  f^(t)  of 
Figure  2.1  and  an  integral  operctor 

Ap(  "  k  ?  [ff(t)  -  f#(t-i)]dt  (2.17) 

where  the  k  is  inserted  to  change  to  metric  units  snd  T  is  the  total 
tracking  interval  (e.g.  12  hr.),  noticing  that  cancellation  takes 


(2.18) 


place, 've  can  re-write  (2.17)  as 


fc(t)dt 


If  is  a  slowly  varying  quantity  (i.e.  it  changes  little  over  a 
period  of  l  seconds) 

m  k4  j^fc(T)  -  fc(0)  J  (2.19) 

12 

Using  the  postulated  drift  of  9  parts  in  10  or  .011  hz  at  S  band 
(2200  Mhz) 

■  80  cm 

Thus  it  is  shown,  that  while  the  high  frequency  variations  from 
the  oscillator  are  much  larger  in  amplitude  than  the  assumed  longer 
terra  drift,  because  of  the  cancellation  afforded  by  the  feed  forward 
loop  of  Figure  2.1,  their  errors  are  much  less  important  than  those 
arising  from  the  slower  (and  invisible!)  drifts  that  take  place  over 
a  longer  tine  scale.  I '.ore  over,  the  effects  of  the  high  frequency  com¬ 
ponents  grow  vithVT;  while  the  low  frequency  variations  go  with  4. 

For  outer  planet  missions,  where  the  light  time  can  be  on  the  order  of 
8  hr.,  the  current  rubidium  oscillators  cannot  perform  adequately 
(they  are  marginal  at  Mars  distances). 


CHAPTER  III 


A  GEOMETRICAL  INTERPRETATION  OP  EARTH-BASED 
TRACKING  DATA  FOR  AN  INTERPLANETARY  SPACECRAFT 


In  this  chapter  our  attention  is  turned  toward  using  the  measure¬ 
ments  of  doppler  and  range  to  determine  the  orbit  of  a  spacecraft. 
These  two  measurement  types,  implemented  as  discussed  in  the  Appendix, 
are  universally  applicable  to  all  phases  of  space-flight  past  the  ini¬ 
tial  powered  flight-Earth  orbit:  trans-lunar,  lunar  orbit  and  lunar 
lander,  interplanetary,  and  soon,  planetary  operations  both  orbitcr 
and  lander.  Consideration  here,  however,  will  be  restricted  to  the 
interplanetary  portion  of  flight. 

This  phase  of  the  mission  is  characterized  by  long  tracking  arcs, 

rather  weak  gravitational  fields  and  hence  seemingly  monotonous  track- 

■  *  1 

lng  geometries.  The  lack  of  rapidly  changing  geometry  such  as  is 
present,  for  example,  in  a  lunar  orbiter  problem,  places  the  most 
stringent,  accuracy  demands  upon  the  tracking  systems  and  physical 
models  sp  that  the  inherent  information  concerning  the  orbit  that  is 
present  can  be  observed  with  the  requisite  fidelity.  Conversely,  this 
same  simple  geometry  can  be  described  by  models  transparent  enough  to 
allow  the  critical  (relevant)  elements  of  the  problem  to  be  discerned. 
The  analysis  to  be  presented  here  has  remilted  from  the  efforts  of 
several  investigators  of  which  this  author  is  only  one.  In  196$, 

Light  (Reference  9)  first  observed  thet  there  was  an  Important 

distinction  between  trans-lunar  and  Interplanetary  flights  due  to  the 

l 

Precedisg  pigs  Mask 
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distances  involved.  In  the  former,  the  angular  diameter  of  the  Earth 
as  seen  from  the  moon  (<^S°)  was  sufficient  so  that  multiple  tracking 
stations,  viewing  the  trajectory  from  slightly  different  aspects 
could  provide  enough  position  parallax  to  fully  determine  the  orbit 
over  very  short  time  intervals.  At  planetary  distances,  this  position 
parallax  is  of  no  consequence,  but  Light  noticed  that  the  rotation  of 
the  tracking  station  about  the  Earth's  spin  axis  impressed  a  diurnal 
signature  on  the  doppler  tracking  data  whose  amplitude  and  phase  were 
of  manifest  importance.  Hamilton  and  Melbourne  (Ref.  2)  formalised 
these  observations  and  put  the  analysis  in  a  form  that  showed  the  dom¬ 
inant  factor  limiting  orbit  determination  accuracy  was  the  accuracy  in 
the  knowledge  of  the  station  locations  themselves.  Later,  this  author 
extended  the  Hamilton/Melbourne  analysis,  which  considered  only  three 
spacecraft  state  parameters,  to  include  the  remaining  variables  and 

9 

produced  a  full-rank  (six  parameter)  description  of  orbit  determina¬ 
tion  accuracy  (Ref.  6).  Perhaps  the  most  extensive  analysis  to  date 
was  published  in  1969  by  this  author  and  8.  R.  Me  Reynold  a  (Ref.  7). 

Zn  this  chapter  we  shall  recall  that  analysis,  its  major  results,  and 
provide  the  background  for  additional  results  in  subsequent  chapters. 

The  Doppler  Data  Equation 

Consider  Figure  3-1  which  shows  the  Earth-relative  geometry  pres¬ 
ent  when  tracking  a  spacecraft  from  a  tracking  station  displaced  from 
the  Earth's  spin  axis  by  the  distance  rg  (it  is  also  elevated  above 
the  Earth's  equatorial  plane,  but  this  is  of  no  consequence).  If 
r»ra  and  the  spacecraft  were  not  moving  with  respect  to  the  Earth, 
the  range-rate  (or  doppler)  would  entirely  be  composed  of  a  harmonic 
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component  of  amplitude 

o)  c°8  6 

where  m  is  the  Earth's  rotation  rate  and  all  other  variables  are  de¬ 
fined  in  Figure  3.1.  If  we  add  the  radial  velocity,  the  first  approxi¬ 
mation  for  the  topocentric  range-rate  is  obtained: 

p  ■  vf  +  o)  rg  cos  6  sin  m  (t-tQ)  ♦  n(t)  (3.1) 

where  n(t)  is  the  data  noise. 

Let  t  be  measured  from  the  time  when  the  right  ascension  of  the  probe 
(a)  end  of  the  station  (£)  on  the  nominal  trajectory  are  in  coinci¬ 
dence,  i.e., 

a  a  _ 
a  -  x  -  0 

In  other  words,  t  »  0  is  when  the  a  priori  estimate  of  the  probe's 
position  will  cross  the  estimated  meridian  of  the  tracking  station's 
coordinates. 

Thus  t  will  define  as  the  time  between  estimated  and  actual 
o 

meridian  crossing.  Evidently, 


where  ■  AX 

ID 

dor  6  o  “  « 

A  A 

AX  X  "  X 

Equation  3*1  becomes 


o  *  vr  +  u>  rg  co®  *  8in  («it  -  a#  +  Ax)  +  n(t) 

■  vy  +  id  rg  cos  *|sin  tut  cos  (Ap  -  t\)  -  cos  „*t  sin 
♦  n(t) 


(3.2) 
(dor  -AX)| 


If  the  a  priori  estimate  of  the  station's  longitude  and  the  probe's 


right  ascension  are  reasonably  accurate,  then 

y*’  y 

f  Ace  -  AXl  «  1 

and  (3.2)  can  veil  be  approximated  by 

p  ■  *fr  *  a  rt  cos  6  iin  «t  ♦  m  r#  cos  ft  (Ax  -  Aa)  cos  ^t  ♦  n(t) 

(3.3) 

Or,  re-writing  in  terms  of  the  unknown  coefficients, 
p  ■  a  *  b  sin  «t  +  c  cos  «t  ♦  n(t) 

Thus,  over  a  time  interval  when  vf.  A,  and  a  can  be  considered  con¬ 
stant,  the  rank  of  the  doppler  information  it  at  most  three.  Equation 
(3.3)  was  derived  in  Reference  2. 

Suppose  we  now  consider  the  three  variables  in  a  first-order  time 
series  expansion,  i.e., 

vr(t)  -  vr  «■  art 
ft(t)  "  t  *  it 
a(t)  »  a  ♦  &t 

Substituting  these  expressions  into  (3>3)  and  using  the  following 
approximation 

cos  (ft  ♦  $t)  m  cos  f  -  sin  ft  (t 

we  obtain: 

p  ■  a  ♦  b  sin  «t  ♦  c  coe  mt  ♦  dt  ♦  et  sin  mi  ♦  ft  cos  mt  ♦  n(t) 

(3.4) 

where 


"a  “ 

“  vr 

b 

m  rt  eog  A 

c 

$ 

to  r,  cos  ft  (ax  -  AOf) 

d 

V 

e 

-to  \  sin  ft  S 

f 

-•«  rg  |  (A\  -  A**)  sin  ft  X  ♦  ^  cos  ftj^ 
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With  the  exception  of  the  d  ten,  all  of  the  coefficient#  are 
written  explicitly  in  terms  of  the  6  unknown  state  parameters;  r,  «, 
a,  vr,  &,  q.  This  d,  or  acceleration,  ten  is  composed  of  two  major 
affects: 

1.  The  component  of  the  8un's  gravitational  acceleration  along 
the  line-of-sight  to  the  probe. 

2.  The  so-colled  "centrifugal  acceleration"  arising  from  a  non- 
to tally  outward  directed  spacecraft  velocity. 

We  shall  leave  the  first  of  these  implicit  for  the  time  being 

but  re-write  d  as 

d  -  a#  (r,  *,  a)  ♦  r(42  cjs2  t>  ♦  42) 

where  a^  is  the  gravitation  component  of  the  acceleration  and  the 

2 

second  tern  is  the  faailiar  m  r  from  elementary  dynaaics. 

Figure  (3.2)  displays  the  doppler  signatures  assuming  the  three 
parameter  model  (a)  and  then  the  six-parameter  model  (b).  In  (b)  the 
signature  from  the  second  day's  tracking  is  shifted  to  the  left  by  2b 
hrs. ,  so  that  the  effects  of  the  presence  of  d,  e,  and  f  can  be 
emphasised. 

The  Ronxe  Data  Equation 

At  t  »  0,  tbs  spin  axis  of  the  Earth,  the  station  vector  r#t,  and 
the  probe  vector,  ?0,  are  co-planar  as  shown  in  Figure  3.3.  Invoking 
the  approximation 

?  ? 

rpl  **  frl 

o(*0)  ■  *(t0)  -  f^tl  coe  (*  -  §>  ♦  n(t0)  (3.5) 

It  can  easily  be  sfcovn  tint  data  t*r<cn  at  other  tines,  tj, 
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Figure  3*2  (a).  Doppler  Signature  Using  First  Three  Terms. 
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DOPPLER  SIGNATURES  ON  TWO  SUCCESSIVE  DAYS  FOR  THE  SITUATION 
SHOWN.  THE  SPACECRAFT  ACCELERATES  OUTWARD.  MOVES  COUNTER 
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DECREASE  IN  THE  DECLINATION  MAGNITlJOE 


Figure  3.2  (b).  Doppler  Signature  Using  Six  Parsaeter  Expansion 


Figure  3.3 


Center  of  Pats  Range  Geometry 


during  the  pass  can  be  referred  to  t  by  adding  to  p(t^)  the  integral 
of  the  measured  floppier  d^t,a  between  points  t^,  tQ  nearly  without  addi 
tional  error.  This  is  because,  over  short  periods  of  time  a  modern 
doppler  system  can  measure  range  change  much  more  accurately  than  can 
the  ranging  system.  A  pass  of  range  data,  then,  (if  it  accompanies 
a  pass  of  doppler  data)  can  be  compressed  to  a  single  point  by  refer¬ 
ring  all  points  to  tQ  and  performing  a  simple  algebraic  average  as  a 
first  filtering  operation. 

The  Concept  on  Data  Compression  as  an  Accuracy  Analysis  Tool 

The  normal  procedure  in  an  estimation  problem  is  to  form  the 
partlals  of  the  data  with  respect  to  the  state  of  the  system  and  es¬ 
timate  that  state  from  the  raw  data  in  a  single  operation.  However, 
note  from  (3.4),  that  the  data  is  written  as  a  linear  function  of  the 
parameters  a  -  f.  Thus,  if  the  system  state  is  first  thought  of  as 
being  these  parameters,  the  accuracy  of  their  estimation  will  not  be 
dependent  upon  the  nominal  vi  lues  since  those  nominal  values  do  not 
appear  in  the  data  partials  (linearity).  This  is  useful.  A  pass  or 
more  of  doppler  data  (perhaps  including  as  many  as  1000  observations) 
can  be  considered  ns  equivalent  to  the  observation  of  the  six  para¬ 
meters,  a  -  f.  Moreover,  since  the  accuracy  of  observation  of  these 
six  parameters  is  not  a  function  of  probe  geometry,  the  computation 
of  how  well  they  cun  be  observed  (for  a  given  data  noise  model)  need 
be  performed  only  once. 

Calculating  the  accuracy  of  the  estimates  of  these  parameters  for 
a  given  length  of  tracking  data  (perhaps  a  twelve  hour  pass)  can  easily 
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be  accomplished  with  any  of  the  estimation  algorithms  detailed  in  the 
Appendix.  Assume  this  has  been  done  and  the  resultant  covariance  of 
errors  is  found  to  be  a  .  Estimation  of  the  state  parameters  r,  6, 

a 

a,  . . it  is  accomplished  by  taking  the  partials  of  a  -  f  with  respect 
to  the  state  parameters  and  using  any  standard  estimation  algorithm. 
Using,  the  normal  equation  formulation,  the  covariance  of  estimation 
errors  is  given  by 


where 


Ar  -  (A^A)*1 


J  $  a"1 
a  Aa 


and  A,  the  partial  matrix,  is  readily  seen  to  be, 


a 

b 

A« 

c 

d 

e 

f 


r 

0 

0 


ft 

0 


-a,  rg  sins 


a 

0 
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-u>  rgcos5 


?SS*2cMW,  £WcoSk.ln?, 

c* 
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?6 

r  Jcosg 

B 

a)  rg  Asinj 


Sa 

0 


vr  * 
1  0 

0  0 


(3.6) 


nr 

0 

0 


2 

0  2r£  2rcos  f(o 


0  -0pr  sins  0 

B 


o)  rsSsins  0  0 


-ci>rgcosf 


(3.7) 


The  gravitational  partials  have  recently  been  put  in  compact 
analytical  form  by  O.ndrasik  (Ref.  8).  They  are: 

^  (2  ’  3  *in2  v) 


. 


u  re  C08  ©[(^j  '  73  )  -  3  77  (r"re  608  x)] 
P  e  P  J 
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■  U  re  008  6  C08  Bin^a  "  a«^[’^r3  "  73)  "  3  “3 
(r  -  re  cos  x)]  (3 


(3.8) 


where 


^  «*  GM  of  the  Sun 

■  Sun  -  probe  distance 
5  ■  Earth-probe-Sun  angle 


r  -  Sun-Earth  distance 

9 


cos  a  ■  cos  cos  jg  cos  (a  -  ag)  +  cos  (,  sin  gg 
X  ■  Sun-Earth-probe  angle 

gg>(yg-  geocentric  declination,  right  ascension  of  the  Sun 

If  it  is  desired  to  study  the  effects  of  augmenting  the  doppler 

data  with  ranging,  the  averaged  effective  p  at  tQ  is  simply  used 

directly  by  augmenting  both  the  Jg  and  A  matrices.  That  is 

r  -2  -1 


Ji  ■ 


(3.9) 


-|?gtjsin(fi-ri  0  0  0  0  0 


(3.10) 


Performing  the  full  analysis  as  indicated  by  (3*6)  will  yield  an 
estimate  of  the  covariance  of  the  state  estimation  errors  for  the  par¬ 
ticular  geometry  used  in  the  evaluation  of  (3.7) .  It  is  useful,  how¬ 
ever,  as  a  first  approximation,  to  inspect  the  behavior  of  Jg^  versus 
tracking  time  and  relate,  op  an  analytic  tttsis,  vhnt  is  now  considered 
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to  be  the  errors  in  the  data  to  errors  in  the  state  estimates. 

Figure  3.4  charts  the  variance  in  the  observation  of  the  parameters 

a  -  f  versus  the  number  of  days  tracking  time  (assuming  a  nominal  12 

hrs./day  tracking  from  a  single  station  and  an  equivalent  uncorrelated 

doppler  data  error  of  1  mm/sec.  for  a  one  minute  sample).  It  is  seen 

that  the  variance  of  the  a  -  c  parameters  goes  roughly  with  l/T,  but 

the  variance  in  d  -  f,  because  of  the  presence  of  the  secular  t  term 

,  .3 

*  in  the  data  equation  goes  with  (l/T)  2. 

Returning  now  to  the  data  partial  matrix,  3.7 ,  it  can  be  observed 
that  three  of  the  state  parameters,  f,,  l  and  &,  are  involved  in  more 
than  one  of  the  data  parameters.  For  short  tracking  spans  (a  few 
days)  and  notably  when  the  nominal  &  and  &  are  small,  the  information 
concerning  each  of  the  state  parameters  comes  largely  from  a  single 
data  parameter.  These  main  partials  are  underlined  in  3.7*  To  the 
extent  that  this  first  approximation  is  correct,  analytic  expressions 
for  the  state  accuracies  as  a  function  of  probe  geometry  can  be  ob¬ 
tained  through  the  information  given  in  Figure  3*4.  These  are 
readily  seen  to  be: 

°r  "  (J3  r3  cos‘T5PS-l.J  +  62cos26  +  $2) 

ofi  -  Ofc/w  !*in  fil 

a  -  ajm  r  cos  ft  (3.11) 

«v  *  °a 
r 

crj  9  rTfl/u  raisin  5 1 

"  ®f/»  rg  co*  6 
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Some  important  geometrical  dependencies  can  be  observed.  Since 
all  spaceflights  to  date  and  those  contemplated  in  the  near  future 
are  restricted  to  lie  nearly  in  the  ecliptic,  the  geocentric  declina¬ 
tion  is  roughly  contained  in  the  interval  -23°  to  +  23°  (from  the 
obliquity  of  the  ecliptic)  and  the  cos  &  term  is  always  near  unity. 
Thus,  to  the  accuracy  of  the  approximation  now  being  used,  the  accur¬ 
acies  of  the  determinations  of  a,  vf,  and  £  are  not  dependent  on  the 
probe  geometry.  However,  lain  5)  can  become  small,  even  zero,  and 
prohibit  the  accurate  determination  of  5  and  &.  This  imposes  an  im¬ 
portant  mission  design  constraint  of  ensuring  that  no  trajectories 
are  selected  which  have  small  geocentric  declinations  during  portions 
of  the  flight  where  rapid  orbit  determination  is  required  (just  prior 
to  planetary  encounter,  after  e  maneuver,  etc.). 

Finally,  in  the  absence  of  ranging  data,  the  r  determination  is 
influenced  by  the  d  partial  as  shown.  The  structure  of  its  magnitude 
is  interesting  and  is  charted  in  Figure  3*5* 

The  following  features  may  be  noted:  l)  The  partial  is  gener- 

J 

ally  large  when  the  Earth-probe-Sun  angle  is  near  l80°  and  has  in¬ 
creasing  intensity  as  the  probe  moves  toward  the  sun;  2)  the  lines 
of  zero  satisfy  the  relationship  that  the  Earth-probe-Sun  angle  equals 
55°  or  125°»  3)  The  partial  is  generally  negligible  for  heliocentric 
ranges  greater  than  1.5  a.u.  and  would  also  be  less  important  for  a 
Mars  transfer  than  for  a  Venus  transfer. 

Of  more  practical  importance  than  lines  where 

£2S.  0 

would  be  oecpciorc  when 
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Figure  3.5.  Mop  of  Pertlal  Derivative  of  the  Differential 
Gravitational  Acceleration  with  Respect  to 
Geocentric  Range. 
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In  this  circumstance,  the  range  could  not  be  determined  by  doppler 
?ac  £ag 

alone.  Since  --*•  2  0  this  can  conceivably  arise  only  when  gr  s  0, 

C-* 

which  from  the  chart  is  generally  in  the  area  between  the  Earth  and 
the  Sun. 

The  Addition  of  Range  Data 

With  these  results  as  background,  it  is  interesting  to  chart  the 
effects  of  augmenting  the  doppler  data  with  range  data.  It  is  clear 
that  when  the  doppler  range  partial  is  weak,  range  data  can  greatly 
benefit  determining  this  component  of  position.  However,  in  favorable 
geometries,  the  range  is  the  best  determined  of  the  position  compon¬ 
ents  with  doppler  data  alone.  The  question  of  whether  a  ranging  sys¬ 
tem  is  warranted  at  all  for  situations  blessed  with  high  gravitational 
fields  certainly  suggests  itself. 

Table  3*1  displays  comparative  uccuracy  estimates  with  and  with¬ 
out  range  for  the  trajectories  experiencing  first  large  (Cose  1)  and 
then  small  (Case  2)  gravitational  range  partials.  In  other  respects 
these  trajectories  are  the  same.  Their  positions  are  plotted  on 
Figure  3.5  as  (l)  end  ©respectively.  When  ranging  data  are  included, 
reliance  can  no  longer  be  made  on  the  approximation  that  a  single 
data  term  determines  a  single  state  parameter.  There  arc  now  impor¬ 
tant  relationships  between  all  the  partials.  The  data  in  Table  3.1 
was  therefore  generated  via  consideration  of  the  full  partial  matrix 
and  algorithms  similar  to  (3.6). 

Several  principles  ore  illustrated  with  these  data: 
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If  * 

I  % 


/■  ■ 

(1)  The  estimate  of  range,  of  cou&e,  imptovea  down  to  the 


••V/J 
-  *•  •  %* 


v  Is  improved 
*v 


the  KPS 


accuracy  of  the  ranging  itself  in  a^^aaes . 

•  **. 

(2)  Even  when  th^‘tcrop'>vsjtocltie$  are  zero,  the  estimate  of 

Jr  *V  *“ 

v  is  improved  by  aiding  rang*  data.  !ftic  improvement  is  slight  for 

n  of  •*'. 

the  KPS  b  53°  case/-,  it  is  $ulte  signijflj&ant  (a  factor  of  2.2  with  two 
days  of  data)  fqfr  the  EPS5-*  l8o  ca&#.*  The  explanation  is  that  when 

/  :l£  *'i;‘ 

the  cross -velocities  artfrlbw  (or<*t&r  )  the  earbh- probe  acceleration 

/  £'  & 

is  caused  dtrictly  by  *Kb  amvitattlonal  field#  and  the  uncertainty  in 


is  caused  ^trictly  by  $fib  |rnvity$|ional  fieldij  and  the  uncertainty  in 

f  •  ^xc,‘  •<  ■ 

the  geoonntaic  rang#feis  tt^pri  ample  reason  /for  a  priori  uncertainty 

j  . 

in  wh#t  the  accef&ration  will  Jmk  Conversely^  supplying  this  range  to 

/  •  m  T 

high  accuracy  #£11  determine  *  tier  geocentric  acceleration.  Because  the 


determine. 


•J  k  a  VH  • 

d  (or  acceleration)  term  if  higt^y  correlate#  with  f  (the  term  that 

•£7%  •  M  t 

.  «  deterr.inetf'v  ).  an  improvement ‘•Jn  knowledge  of  d  brings  with  it-  an  im- 
f  a  •  •  j,  f 

£  proven^**  1  ft  f,  and  the  agcuracj^>f  1m  consequently  improved.  In 

*■  fact,  -.if*  d  Js  perfectly  uncertainly  in  f  will  improve  by  a 

ji?  •*...  •  •  J  €  : 

faetj^c  ofj2.P,  tH@  same  Sdprov^pnt  already ^observed  in  v^  for  the 

£ca^jc.  The  full  2^Yimpf^bment  4s  nol^expcrienced  for  EPS  «  53° 

se,  Jn  the*vfcinit#  of  1J&  sun  (<  1.5  W),  a  low  fca/gr  really 

Ay'  .  •  #»  •  v 

•bft  .  •  •  •  •  **  Sgi 

.-»>  inqpj.es  qpt  that  or  is  .ErnallSlut  that  it"id,perpendicular  to  the 

*•  *  *  ;r  r 

JLine  oS*sight.  Titos  unpertaii^es  in  the  'Jjross  positions  produce 

^gnifianp t  variations  d£  ranging  $^ply  will  not  remove  this 

.  ,  uncertainty.  •  f  Jr|?  •  1 

(3^*  No  correspond jhg  l^rovement  ocqui-s  in  v  when  the  a  prior 


&  - 

•*V- 

.  —  >*. 


V  (3i.*Ko  correspond *ng  iojA-ovemcnt  occ^i-s  in  v  when  the  a  prior 

m  %  ,7*  .||  .J  6 

v^"^|^zero.*‘ ^ftese  <Hb*e  te^pVhich  detenfljines  v  is  uncorrelated  with 
A  - *•  - - -  aijpply. 


d  and 


ced  does  not 


•  1 

(U)  When  the  cross-velocity  is  non-zero,  the  addition  of  range 


data  significantly  improves,  the  determination  of  this  velocity.  The 
reason  here  is  clear.  Supplying  range  data  makes  the  entire  power  of 
the  d  term  available  to  improve  the  velocity  estimation.  This  improve* 
ment  occurs  only  in  a  single  direction  -  along  the  a  priori  cross- 
velocity  itself;  c.g. ,  notice  that  for  the  cose  where  both  v  and  v 

6 

are  3  km/s,  little  improvement  occurs  in  either  individual  direction  - 
the  improvement  is  experienced  in  the  direction  45°  from  the  axes  we 
have  chosen.  Coupled  together,  Items  (2)  and  (4)  explain  the  often 
observed  but  previously  curious  result  that  the  addition  of  range 
data  considerably  increases  the  rapidity  of  orbit  determination  follow¬ 
ing  a  midcourse  correction  where  only  the  velocities  have  been 
corrupted. 

(5)  Similar  computations  were  made  to  demonstrate  the  dependence 
of  1  -  4  on  the  precision  of  the  range  data.  The  velocity  improvement 
was  found  to  be  unrelated,  however,  as  long  as  the  measurement  was 
considerably  (X10)  more  precise  than  the  accuracy  of  range  determined 
from  Doppler  alone. 

In  Summary,  range  data  improves  the  accuracy  of  a  single  com¬ 
ponent  of  position  as  expected  and,  in  addition,  one  component  of 
cross-velocity.  The  former  is  strictly  a  function  of  the  measurement 
itself,  the  latter  depends  upon  the  magnitude  and  direction  of  the 
cross-velocity. 

1  ,,  *  t  - 

v  . 

o 
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Table  3.1.  The  Influence  of  Ranee  Data  on  Orbit  Accuracy 


2  daya  of  data 


CASES 

or.  *■  «/•« 

0n>,  */aec 

EPS  -  l80°  (2ag/>r  -  lU  x  lO*1*1  see'2) 

r 

t* 

■ 

0 

lb  ranging 

9.9  0.13 

0.23 

■ 

0 

Ranging 

0.01  0.058 

0.23 

Ta 

• 

3 

lb  ranging 

193.  0.11 

0.23 

© 

r6 

m 

0 

Ranging 

0.01  0.0027 

0.23 

r 

or 

m 

3 

Ho  ranging 

373.  0.11 

0.19 

m 

3 

Rangier 

0.01  0.11 

0.11 

EPS  »  53°  (»ag/3r  -  0.25  *  10*lU  aec'2) 

▼ 

or 

m 

0 

Ho  Ranging 

389.  0.13 

0.22 

T. 

- 

0 

Ranging 

0.01  0.096 

0.21 

© 

V 

or 

• 

3 

Ho  ranging 

IU13.  O.lU 

0.2 2 

■ 

0 

Ranging 

0.01  0.0039 

0.22 

Trajectory  charactcriatics:  r  ■  26  r.  lO^kaj  ■  0; 

>  6  •  15^5  ft;  * 

0 

1  a/iec  o  reneing  -  10  ti;  saople  rate  ■  l/tvo  daya; 

aacple  rate  ■ 

l/aln* 


TU  EFFECT  OF  KAIIOi  I/XATIOff 


i>  *<  ■  •  >  y i  ■»  y  (•’  W  ' I  t.th 


In  this  chapter  we  provide  a  description  of  the  errors  Introduced 
to  the  state  estlaates  bp  the  preaencc  of  1)  errors  in  specifying  the 
poaltlona  of  the  Earth-based  tracking  stations  and  2)  snail  forces 
acting  on  the  spacecraft  and  perturbing  Its  notion.  The  data  con- 
press  ion  point  of  view  developed  In  the  previous  chapter  is  used  and 
both  tbs  situations  where  dopyler  data  alone  are  enployed  and  then 
when  ranging  data  aegnent  th*  basic  doppler  observable  are  consider* 
ed.  The  analysis  presented  here  will  be  accurate  for  only  abort  arcs 
cf  tracking  data  (perhaps  one  to  two  weeks)  but  will  provide  insight 
to  that  important  case  and,  in  addition,  assist  in  the  interpretation 
of  the  long-arc  analysis  presented  In  Chapter  VI. 

Both  the  station  location  errors  and  unnodclcd  forces  enter  the 
doppler  observable  equation  in  relatively  sinplc  ways.  Reconsider 
Equation  3.b.  The  pertlals  of  the  conpressed  data  (a  through  f)  with 
respect  to  the  errors  to  be  considered  are: 


Preen*  pip  im 


17 


where  r  A  the  tracking  station'*  distance  froa  tbs  Barth's  spin 
*  axis 

X  ft  the  longitude  of  tbs  tracking  station 

*n«  4  «.  KCImtU.  1.  U»  Mm*.  41m.!. 

Rots  that  tbs  thi*d  component  of  tbs  station  location,  the  distance 
above  the  Barth's  equatorial  plane,  and  the  reasoning  two  conponents 
of  acceleration  orthogonal  to  the  line  of  sight  do  not  enter  this 
analysis. 

The  p  rob  lea  becomes*  knowing  the  effect  of  the  aodel  errors  on 
the  data  observable,  coapute  the  relationship  between  data  errors  snd 

state  estlaates. 

Inversion  of  the  Doppler  Partial  Derivative  Matrix 

In  the  previous  chapter,  no  explicit  advantage  was  taken  of  the 
fact  that  the  original  data  had  been  compressed  into  an  observable 
vector  whose  d  lac  ns  ion  equaled  that  of  the  state  vector.  That  is,  if 
t  a  ■  A  i  x  ♦  n 

Instead  of  forming  the  state  estimate  by 

I  ft  •  (A1  W  A)’1  AT  If  g  a  (h.2) 

we  can  slaply  write 

ft  x  •  A"1  6  a  (h.3) 

since  A  is  an  invertible  aatrix  and  is  given  by  Equation  3.7. 

The  simple  structure  of  this  aatrix  allows  its  analytic  inverse 
to  be  taken  with  relative  ease 


ha 


•In*  co« g  ♦  «-  -  2rl  gj  ♦  2rv$  *ii*§  coag  «rr  2rJ  2r  coag6 


where  the  partial!  *7,  are  given  by  (3.8). 


Rotation  to  Cartesian  Coordinate! 

Pre-multiplication  of  (4.1)  by  (4.4)  will  yield  the  desired  re* 
suit  of  the  sensitivity  of  the  state  estimates  to  variations  in  the 
three  model  parameters  under  consideration.  However,  in  order  to  in¬ 
terpret  the  results  and  propogstc  forward  in  an  inertial  frame,  a  ro¬ 
tation  from  the  actual  trajectory  related  system  being  used  to  an 
inertial  cartesian  systca  is  required.  This  is  easily  shown  to  be: 


Where,  for  example,  r  is  the  position  deviation  in  the  direction  of 
increasing  declination  measured  from  the  nominal  trajectory.  All 
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quantities  on  the  left  are  in  an  inertial  frane,  whereas  quantities 

on  the  right  are  measured  from  the  actual  probe  state.  For  example, 

gv  ia  the  error  in  the  probe's  radial  velocity;  gv  is  the  error  in 
rT  rI 

the  velocity  coxponent  along  the  nosiinal  radial  direction,  which  is 

lnertlally  fixed. 


Analytic  Sensitivities  to  Model  Parameters 


The  doppler  analysis  can  now  be  completed.  First,  pre  Multiply¬ 
ing  (4.1)  by  (4.4) 


U-m 

*>r. 


.2 


2r _ lg  sin&  cosg 


r  ^ 
»jr 


WjlV 

♦  arV2rr 


r 


tang 


-L 


tang  r# 

0 

0 


Vr,(tan"2g  ♦  1) 


2  n 


«/r. 


gd/?r 


0 

0 


«*»»  -  i*  -  i*  co.2 


(4.6) 


(4.7) 
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and 


0 

0 

0 

0 

0 


<4.8) 


Tha  alieht  transformation  fro*  d\  to  r#d*  haa  been  Mde  ao  that 
both  components  of  a tat ion  location  error  can  be  expressed  In  ka. 

Rote  from  (4.6)  and  (4.7)  that  the  first  immediately  apparent  effect 
is  that  an  r#  ei-ror  produces  a  perturbation  in  declination,  while  a 
longitude  error  produces  a  right  ascension  perturbation  in  the  est las¬ 
ted  state.  These  effects  were  first  demonstrated  by  Hamilton  and 
Melbourne  (Reference  2 ) .  In  addition,  an  r#  error  Introduces  a  range 
error  whose  magnitude  is  dote  rained  by  a  rather  complex  interrelation¬ 
ship  between  elements  of  the  nominal  state.  This  error  can  be  ex¬ 
cessively  large  in  certain  circles tnnees,  particularly  when  gd/jr 
becomes  small  or  vanishes.  Also,  there  are  small  velocity  perturba¬ 
tions  as  well. 

The  state  errors  introduced  by  station  longitude  errors  also  in¬ 
clude  a  range  term  which  depends  strictly  upon  the  ratio  of  the  geo¬ 
centric  acceleration  parti  a  Is  with  range  to  those  of  right  ascension. 
Finally,  an  unmodeled  acceleration  error  produces  sn  error  in  range 
only  as  shown.  Note  that  all  three  error  sources  produce  large  range 
errors  when  jd/jr  is  small.  Singularities  occur  in  several  state  com¬ 
ponents  for  both  r#  and  \  errors  when  the  probe  declination  is  aero. 
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In  rotating  to  the  inertial  cartesian  coordinate  system,  one 
further  substitution  must  be  made  to  keep  the  expressions  obtained 


from  getting  unwieldy.  Define 

a  (iA.  i*/ 

k  A  1  *inft  co*6  4  ar  \2r/ 

tang 


-  (J2  ♦  CO.2,) 


Then  using  the  transformation,  (4.5),  we  obtain  the  sensitivities  in 
the  plane-of-sky,  inertial,  cartesian  frame. 


2r/_  2dk 

»  #f 


S-4*  r  w. 


t»n6  rf 


“ rJ/tsnjr # 

2rAk  jr _ 

r,  ad/a*’  *  tangr# 


(4.9) 


+  ^  (tan'2g  «■  1) 


Procosfk  rising  2rc  on  ter 

r,  8d/dr  "  tangr#  +  r# 

The  longitude  sensitivities  sre  much  simpler: 


-rcos 


(4.10) 


(l± /l±  r  h 
/ai  /a*  \ 

(6  cosf )  /  *r  rfly 


v,cosg 


l/(&d/&r) 


0 

0 

0 

i/ad/t'r 

cr  cos^/?d/ar 


Together  (4.6)  through  (4.11)  provide  o  rather  complete  descrip* 
ticn  of  the  effects  of  these  model  errors  over  short  tracking  periods 
with  doppler  only  data  in  both  coordinate  frames  of  major  interest. 
NUserical  examples  are  deferred  until  ranging  data  are  re-introduced 
to  vhe  problem. 

The  Effect  of  Ranging  Data 

When  range  data  are  employed  the  problem  can  no  longer  be  worked 
in  the  simple  analytic  fashion  „ust  demonstrated.  Wc  begin  by  pz-c- 
senting  a  simple  approach,  done  with  the  square-root  estimation  al¬ 
gorithms  discussed  in  the  Appendix,  for  the  computation  of  the  sensi¬ 
tivity  of  the  least-squares  estimator  to  the  modeling  errors  under 
consideration.  In  general  notation,  the  data  set  ge,  is  related 
linearly  to  the  state  through  the  familiar  data  equation, 
ga  ■  A  ax  ♦  n 

where  n  is  the  additive  noise.  In  this  case,  the  A  matrix  is  given 
by  (3*10)  and  the  inverse  covariance  of  the  noise,  J,  is  given  by 
(3*9)*  Thus  pre-nultiplying  through  by 


An  orthonormal  transformation  P  that  places  J*A  in  upper  triangu¬ 
lar  form  is  found  and  the  indicated  estimation  procedure  as  outlined 
in  the  Appendix  is  executed 


Thus  the  upper  n  x  m  portion  (n,  m  is  the  dimension  of  the  state,  data 
vector;  in  this  case  6,  7)  of  PJ^  specifies  the  relationship  between 
data  noise  and  estimation  error.  Finally,  if  n  is  replaced  by  the 
m  x  3  matrix,  dn/fcr  ,  r  ,  a  ,  the  sensitivity  matrix  of  errors  with 

•  V 

respect  to  the  three  parameters  being  studied  -  station  locations  and 
non-gravitational  forces  -  is  found. 


The  bottom  row  gives  the  change  in  the  residuals,  i.e.,  the  component 
of  the  data  orthogonal  to  the  column  space  of  J^A,  with  the  parameter 
variations.  This  is  a  single  number  per  parameter  and  conveniently 
gives  the  norm  of  the  total  residual  vector  caused  by  a  unit  parameter 
perturbation.  This  lost  comment  is  made  parenthetically  but  we  shall 
return  in  Chapter  VII  to  what  is  now  a  vague  point  to  formally  gener¬ 
ate  machinery  that  will  allow  computation  of  the  increase  of  the 
weighted  sum  of  squares  of  the  residuals  caused  by  the  presence  of 
unsodeled  parameters. 


Numerical  Example  -  Response  to  Unmodeled  Forcci 

In  order  to  apply  the  foregoing  analyaia  to  a  realistic  situation 

a  set  of  Viking  '75  trajectories  is  chosen.  There  are  13  cases  in  all 

« 

spanning  the  launch  -  arrival  date  space  currently  planned  for  the 
mission.  All  but  4  of  these  cases  uce  nominal  initial  conditions  at 
30  days  from  Mars  encounter.  On  two  separate  trajectories,  conditions 
at  encounter  minus  10  days  and  encounter  minus  60  days  are  also  ana¬ 
lysed.  The  reason  for  the  multiplicity  of  cases  is  to  attempt  to  de¬ 
termine  if  there  are  large  variations  in  the  answers  obtained  over 
the  geometry  excursions  typical  of  a  single  mission.  We  wish  to 
study  the  influence  of  arcs  considering  l)  doppler  only  data  and 
2)  doppler  plus  range  data.  Table  4.1  identifies  the  case  numbers 
and  lists  what  will  prove  to  be  the  single  most  important  variable, 
geocentric  declination. 

Table  4.2  summarises  the  pertinent  results  from  this  study.  In 
accordance  with  Equation  4.8,  the  doppler  -  only  results  show  pertur- 
bation  in  the  r  coordinate  only.  For  a  10"  km/sec  applied  force, 
the  results  are  as  shown  in  the  first  column.  There  is  little  vari¬ 
ation  over  the  trajectories  inspected  from  a  low  of  435  km  to  a  high 
of  5**5  km.  The  reason  is  clear:  the  error  depends  largely  on  the 

22a 

term  which  is  determined  by  the  solar  gravitation.  Since  Mars  is  at 
a  roughly  fixed  distance  from  the  Sun,  the  only  variable  is  the  EPS 
angle  (Equation  3.8)  and  this  does  not  vary  greatly  over  the  arrival 
window. 
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Tabic  4.1 

Vikang  '75  Trajectories  Selected  for  Smell-Force  Study 


Case  No. 

Launch  Date-Arrival  Date 

Time 

6  (<*««.) 

I 

September  9,  '75 

E  -  30  days 

11.2 

August 

8,  *76 

II 

August 

14,  *75 

E  -  30  days 

14.5 

July 

24,  *76 

III 

Aiigust 

19,  ’75 

E  -  10  days 

.48 

September 

2,  '76 

IV 

August 

19,  '75 

E  -  30  days 

5.1 

September 

2,  '76 

V 

Augubt 

19,  '75 

E  -  60  days 

11.4 

September 

2,  '76 

VI 

September 

3,  '75 

E  -  30  days 

8.8 

August 

11,  '76 

VII 

September  13,  '75 

E  -  30  days 

11.9 

August 

2,  '76 

VIII 

August 

16,  '75 

E  -  30  days 

18.5 

July 

1,  '76 

IX 

August 

June 

1,  '75 

6,  '76 

E  -  10  days 

19.9 

X 

August 

1,  '75 

E  -  30  days 

22.0 

June 

6,  '76 

XI 

August 

1,  '75 

E  -  40  deys 

23.6 

,  i 

June 

6,  '76 

XII 

August 

August 

14,  '75 

13,  '76 

E  -  30  dnys 

10.1 

12.4 


XIII 


August  14,  '75 
August  3,  '76 


E  -  30  days 


Table  4.2 


Results  for  Viking  '75  Small  Force  Study 


_  11  P 

Assumed  Force  of  10  km/sec  in  Radial  Direction 


Case  11 

Doppler 

Only 

Doppler  $ 
Range  11 
Days  of 
Tracking 

U  (tan) 

ARA  km 

th  m/sec 

m/sec 

Residual 

6 

1 

538 

21.3 

(18.) 

-13.4 

(-11.) 

-.046 

(-.038) 

.031 

(.025) 

.436 

11.2° 

II 

545 

15.2 

(12.3) 

-16 

(-13) 

-.032 

.036 

.48 

III 

515 

23 

(23) 

8.6E-3 

(8.6E-3) 

-.1 

(-.1) 

P ,  .  * 

.  IE-1 

■ 

IV 

535 

36.8 

(35.9) 

ESI 

-.082 

(-.075) 

.284 

a 

V 

555 

21.6 

(18.1) 

-13.6 

(-11.4) 

-.047 

(-.034) 

mmwm 

.44 

mi 

VI 

545 

27.7 

(24) 

-10.5 
(-  9.2) 

-.058 

(-.052) 

WUk 

.386 

■ 

VII 

545 

19.3 

(16.0) 

bum 

-.042 

(-.035) 

.032 

(.027) 

.45 

11.9° 

VIII 

535 

8.9 

(6.99) 

-18.6 

(-14.6) 

-.0195 

(-.0:5) 

.044  . 

(.0346) 

.52 

SB 

IX 

540 

7.7 

(5.99) 

-19.3 

(-15.) 

BH 

.0**4 

(.0344) 

.537 

19.9° 

X 

520 

mm 

-20.7 

(-15.6) 

-.01 

(-.0078) 

.04? 

(.037) 

.568 

22° 

XI 

435 

1.34 
(  .94) 

-22.8 

(-16.1) 

-.002 

(-.0014) 

mm 

.648 

23.6° 

XII 

543 

24.4 

(20.8) 

-12.0 

(-10.4) 

-.052 

(-.045) 

.027 

(.024) 

.414 

10.l£ 

XIII 

m 

19.3 

(16.0) 

woo 
•  • 

H 

H  H 
l  • 

-.041 

(-.034) 

1  ...  _  » 

.032 

(.026) 

.455 

■ 
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Insertion  of  ranging  data  dramatically  changes  the  situation. 

Note  that  the  *>00  km  perturbation  in  position  is  now  reduced  to  appro¬ 
ximately  30-40  km  distributed  between  the  right  ascension  and  declina¬ 
tion  directions.  In  addition,  velocity  perturbations  arise  in  the 
a,  $  directions  as  shown.  These  later  parameters  are  variable  but  the 
maximum  obtained  is  .1  m/sec  in  the  $  direction  for  case  III. 

Thus,  if  position  determination  is  the  main  criterion  used  for 
evaluating  the  orbit  determination  accuracy,  adding  the  raige  point 
greatly  benefits  the  accuracy  in  the  vicinity  of  the  data  being  used. 
The  orbit  will  not,  however,  propagate  as  well  because  of  the  appear¬ 
ance  of  velocity  errors.  If  a  gravitational  field  free  propagation 
model  is  used  as  a  rough  indicator  of  how  these  errors  will  propagate, 
the  .1  m/scc  error  occurlng  in  case  III  will  propagate  to  500  km  in 
approximately  50  days.  Of  course,  it  is  in  a  different  direction  than 
the  original  500  km  error,  which  was  along  the  Earth  line  of  night. 

The  numbers  in  the  parenthesis  are  the  perturbations  arising  when 
the  acceleration  is  included  in  the  estimate.  To  obtain  these  numbers, 
an  a  priori  uncertainty  of  10  km/sec  was  Included  in  the  model. 

The  ability  to  solve  for  the  acceleration  and  decrease  the  resultant 
error  varies  from  case  to  case,  but  is  nowhere  particularly  powerful. 
The  residual  column  displays  the  r.s.s.  of  the  weighted  residuals 
after  the  fit  caused  by  the  unmodelcd  acceleration.  Since  the  ex¬ 
pected  sum  of  squares  from  data  noise  alone  would  be  unity  (this  will 
be  shown  in  Chapter  VII),  this  is  not  a  dramatic  Increase  in  any  of 
the  cases. 

The  doppler  plus  rnn^e  results  are  P  function  of  all  the  nominal 
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p.obe  coordinates,  but  a  glance  at  Table  4.2  lndlcatea  a  atrong  depen¬ 
dence  on  geocentric  declination.  To  aee  this  more  clearly,  Figure  4.1 
plots  the  aagnltudo  (note  that  the  sign  between  each  position  and  its 
corresponding  velocity  coordinate  la  reversed  in  all  cases)  of  the 
perturbations  versus  this  factor.  The  dependence  is  narked  and  quali¬ 
tatively  describes  as:  for  low  declination,  the  perturbations  enter 
the  g  and  £  directions  -  as  declination  increases  this  is  traded  for 
errors  in  the  Q  and  ^  directions.  The  reason  for  this  behavior  is 
clear.  Since 


and 


& 

36 


sing 


as  g  0,  g  and  &  ccn  be  perturbed  to  reaovc  the  d  tern  residual  with¬ 
out  introducing  large  changes  in  the  b  and  e  residuals.  This  is  borne 
out  by  the  residual  column  of  Table  4.2.  The  topocentrlc  terns  involv¬ 
ing  a  and  ^  contain  no  such  geometric  dependence  and  are  relatively 
constant  froa  case  to  case. 

In  susnary,  the  main  conclusions  reached  here  for  the  short  arc 

are: 

1.  Small  forces  can  induce  relatively  large  errors  in  range 
when  using  doppler  only  data. 

2.  Ranging  data  removes  this  error  nearly  completely,  but  at 
the  expense  of  errors  in  g,  0,  ^  & 

3*  The  position  determination  is  much  improved  in  the  vicinity 
of  the  data  taken  when  this  includes  range,  but  the  orbit 
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Figure  fc.l.  Estlcction  Errors  Introduced  froo  s  10"11  kn/sec2 
Acceleration. 
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?7  ONV  JV 


^ined  d<*»  i'0*'  Propagate  „g  veil.  For  the  Viking  cases 
•^led,  tl*  brc®k*«ver  point  is  on  the  order  of  50  days 
ProPap>/lti°n. 

****  kbil***  to  e°lv,s  'or  the  force  over  a  short  tisw  period 
i*  *odest* 


A  SQUARE- HOOT  OtOOAUZED  OOHSnt*  0PT20N 


The  p  reeved  log  chapter  explored  the  defredatloe  la  filter  per* 
fonaence  due  to  the  existence  of  tvo  limiting  error  ao ereeet  1)  sta¬ 
tion  location  errors,  and  2)  —11  accelerations  perturb  lag  the  notion 
of  the  probe.  Both  developments  were  only  valid  for  relatively  short 
area  of  tracking  data.  Za  this  chapter,  we  wish  to  develop  two  tech¬ 
niques  for  calculating  filter  perforaance  la  laperfsetly  node  led  situ¬ 
ations.  They  will  be  presented  here  as  self-contained  studies  and  as 
apparent  diversions  fron  the  aala  pursuits  of  the  previous  chapter, 
but  will  both  1  mediately  be  applied  to  the  problem  of  filter  response 
to  poorly  modeled  environments :  This  tine  over  long  tracking  data 
spans.  Ve  wish  to  develop  these  techniques  as  an  application  of  the 
Householder  transformation  approach  to  the  solution  of  the  linear 
least  squares  problem.  This  technique,  which  is  sometimes  referred  to 
as  square  root  filtering,  originally  came  from  a  method  by  Householder 
for  inverting  square  matrices  (by  rotating  them  to  upper  triangular 
form  and  Inverting  with  backwards  substitution  -  Reference 9).  This 
was  extended  to  non-square  a •  trices  by  Businger  and  Golub  (Reference 
10).  Later,  Hanson  and  Lawson  identified  two  basic  algorithms  to  1) 
construct  a  Householder  transformation  matrix  in  compactly  stored  form 
and  2)  to  multiply  a  vector  by  such  a  matrix  (Reference  11).  Dyer  and 
McReynolds  derived  sequential  formulae  for  adapting  the  square  root 
concepts  to  Knlman-Bucy  Altering  Including  the  optimal  treatment  of 


process  noise  (Reference  1?). 


A  short  derivation  of  the  use  of  these  techniques  for  the  slm- 
plest  of  estimation  problems  is  presented  in  Appendix  HI  where  it  say 
be  cor  pared  to  the  classical  least  squares  and  sequential  formulae. 

The  foil  owl  nr  assumes  faallisrity  with  that  material.  We  do  not  deal 
explicitly  vith  the  numerical  advantage  accrued  by  using  the  square- 
root  point  of  riev;  the  material  in  References  9*12  is  thought  to  be 
conclusive  on  this  point.  Rather,  it  is  assumed  that  any  filtering  or 
evaluation  problem  which  can  be  east  in  this  framework  benefits  ade¬ 
quately  to  Justify  the  effort. 


the  Square  Root  Consider  Option 


Dyer  and  McReynolds  (Reference  13)  have  extended  the  application 
of  the  Householder  algorithm  to  include  the  case  where  some  of  the 
unknown  parameters  arc  not  included  in  the  estimation  set.  It  is  de¬ 
sired  to  calculate  the  contribution  of  these  parameters,  referred  to 
as  "consider"  parameters,  to  the  resultant  description  of  the  estima¬ 
tion  errors.  We  recall  chat  development  briefly  here.  Suppose  that 
an  information  array  has  been  obtained  in  terms  of  two  classes  of  para¬ 
meters,  x  and  y.  That  is 
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Suppose  that  Instead  of  estimating  x  and  y,  x  only  is  estimated  end  y 
is  left  at  the  a  priori  veluc  of  zero.  Setting  y  to  zero  and 


pre-multlplylng  (5.1)  by 
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Transposing  the  x  in  (5*2)  to  the  left  bend  side,  poet -Multiplying  the 
resultant  equation  by  its  transpose  and  taking  ensemble  averages  yields: 

E[<*  -  ,)(4  -  x)T]  t  ^  .  R,- V1**  (5-3) 

where 

/,  « t{r  /] 

That  is,  y  is  thought  of  as  a  random  variable  (uncorrelated  with  x) 
with  covariance  Whether  the  resultant  ea  time  tor  of  x  is  unbiased 
will  depend  upon  the  mean  of  y:  this  is  usually  assumed  as  zero. 
Equation  (5*3)  can  easily  be  shown  to  be  equivalent  to  the  classical 
least  squares  result  for  the  same  circumstances  as  derived  in  the 
Appendix. 

Generali  ring  to  Incorrect  Anri  or!  Information  and  Cross  Correlation 
Between  the  Kstlr.oted  and  Considered  Parameters 

Zn  both  the  simple  circumstances  of  estimating  all  the  relevant 

variables  and  considering  a  subset  of  those,  it  has  been  assumed  that 

the  opriori  covariance  used  in  the  filter  design  for  the  estimated 

variables  hr.s  been  correct.  In  addition,  t£e  r.priori  cross 
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correlation  between  the  eatlxnted  and  consider  variables  baa  been 
assumed  tero.  Other  lnvestlfAto,-s  have  studied  the  effects  of  in¬ 
correct  a  priori  on  the  esti  fitted  parameters  both  from  the  batch  aid 
sequential  filter  points  of  view  (See  Song,  Referenced,  and  Rishimura, 
Reference  15). 

It  can  be  useftil  in  several  ways  to  combine  the  notions  of  in¬ 
correct  a  priori  vith  the  consider  option  to  aid  in  practical  filter 
design.  As  a  concrete  example,  suppose  thot  two  parameters  exist  that 
enter  the  data  in  nearly  a  linearly  dependent  manner.  It  would  seem 
reasonable  as  a  practical  matter  to  estimate  only  one  of  these  vari¬ 
ables,  but  input  to  the  filter  an  a  priori  uncertainty  high  enough  to 
(partly)  corpensatc  for  not  estimating  the  second  variable.  In  order 
to  accurately  evaluate  the  effects  of  doing  this,  we  need  to  reflect 
that  ve  have  used  incorrect  (too  high)  a  priori  uncertainty  on  the 
first  variable  and  have  not  estimated  the  second.  This,  and  other 
similar  notions,  can  be  treated  more  rigorously  if  any  apriori  cross 
correlation  between  the  x  and  y  variables  can  be  handled  by  the  de¬ 
rived  algorithms. 

The  exact  problem  to  be  investigated  is  the  following:  the 
filter  will  estimate  the  x  variables  and  will  include  apriori  data 
regarding  these  variables  whose  uncertainty  is  specified  as 

*{<4  -  x)(i  ■  X)T].  priori  •  \ 

The  actual  data  is  influenced  by  on  additional  set  of  variables,  y. 
Moreover,  a  i*  not  necessarily  an  accurate  model  of  the  apriori 
data.  Hie  true  situation  is 


-  *)(*  -  .)*]  Hi  -  x)(y  -  y)TJ 
V  -  y)(ft  -  x)T]  l£(y  -  y)(>  -  y)1] 


4r  N 


i  priori 

Additional  data,  relating  to  both  x  and  y  la  takan  and  adjoined  to  the 
aaaured  a  priori  data. 
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where 


ia  tha  aprlori  aatlaaU  of  x  having  d  Irena  Ion  n 
a  ia  tha  ■  d Irena  local  data  vac  tor 

y  la  tha  p  direnaional  act  of  cocalder  variablca 

*ap,a  ia  tha  arror  in  and  a  respectively 
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•M 

Tha  total  data  vector,  a^,  la  weighted  by  the  inverse  square -root 
of  tha  assured  arror  covariance, 


yielding, 
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Equation  (5.9)  la  triangular!***  by  an  (n«u)  a  (nna)  orthonorual 


■atrlx,  P,  resulting  In 


:  a  *  . 
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(5.6) 


ubarr  P  haa  been  partition**  as 

n  ■ 


Recognizing  that  the  estimate  of  x  is  alaply 

*  „-l  " 

x  -  Rx  zx  , 


(x  -  x)  «  RxX  {r^  ♦  Pu^eap  +  PigA,*  «} 


(5.7) 


Multiplying  (9.7)  by  its  transpose  and  taking  ensemble  averages  results 
in  the  posteriori  covariance  on  x, 

/(S  -  *><*  ■*?*£•  ♦  v>  ‘Ip]\fpn 

+  Pll\t  <«.ppT]  Rxy  *  pll*^  *Ce«P*»p]  ^  ’ll 
♦  I e,1]  a;*tp^  j  (5.8) 


Equation  (5.8)  must  be  Interpreted  in  terras  of  the  original  problem 
statement  to  arrive  at  the  final  expression.  Since  y  is  zero, 
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Invoking  •  consequence  of  the  orthonormolity  of  P,  can  be  written 


in  tern  of  P^  as , 

V*  ”  1 


p  p  ^ 

*11*11 


Substituting  the  above  in  (5.6)  yields  the  desired  result, 

>T 


-  R 


'xy^xy^x*  '11 


<  •  { "x/yKxy 

-  Pn\|  fyly  +  Pl).  r'x*  *x  \f  '  J]  Pll  4  1  }  "i1*  <5-9> 

Finally,  the  posteriori  eovarinnee  between  x  and  y  is  readily  seen  to 
be: 

Ef(i  - ,)(}  -  y)1]  «  -  "x1  V  V  pu  \J  <5-l0> 

The  P  matrix  is  not  normally  explicitly  calculated  in  the  stand¬ 
ard  Householder  transformations  and  special  provisions  must  be  made  to 
calculate  P^.  R  Hanson  has  poink-ed  out  in  a  private  communication 
that  the  most  convenient  way  of  accomplishing  tills  is  to  adjoin  an 

(n<m)x  n  matrix  to  the  weighted  data  vector,  ,  of  the  following 

*t 

form: 

P 


‘i  ill  r *i'|  pu ' 

M  6  J  [vp;r 


P^  (as  well  as  will  appear  as  a  result  of  this  operation.  It 
does,  however,  requi  re  new  storage  requirements  as  shown,.'  although  if 
that  causes  a  problem,  can  be  fonned  sequentially  and  Pg^  never 
need  be  explicitly  stored. 
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A  Stochastic  Square-Root  Option 


Perhaps  the  most  important  advantage  of  sequential  filtering  over 
the  classical  batch  techniques  is  its  ability  to  deal  optimally  with 
stochastic  forcing  functions  in  the  state  dynamics.  The  simplest  of 
circumstances  would  be  to  have  a  state  equation  of  the  form 
x  «»  Fx  +  v 


or 


i+1 


Fx,  +  w. 


in  a  discrete-time  system,  where  w  (w^)  is  a  continuous  (discrete) 
white  (sequence)  process.  In  this  circumstance,  the  filter  does  not 
explicitly  estimate  w  but  uses  the  covariance  of  w  in  the  matrix 
Ricatti  equation  for  the  covariance  of  the  state  estimation  errors  to 
make  a  continuous  calculation  of  the  optimum  Kalman  gains  possible. 

In  the  case  of  correlated  process  noise,  en  explicit  estimate  of  w 
may  be  made  by  augmenting  the  state  equation  to  include  the  time  cor¬ 
related  noise.  The  differential  equation  governing  the  correlated 
noise  will  be  itself  driven  by  white  noise  (i.e.,  we  ore  dealing  with 
a  linear  shaping  filter). 

It  is  useful  to  calculate  the  effects  on  filter  performance  due 
to  unmodeled  stochastic  processes.  That  is,  the  filter  is  ignorant  of 
the  existence  of  these  processes  and  i3  optimal  with  respect  to  a  re¬ 
stricted  model  only.  However,  the  true  filter  performance  in  the 
presence  of  these  stochastic  processes  mt.st  be  calculated. 

This  procedure  might  be  termed  a  stochastic  consider  option.  In 
general,  the  process  will  affect  both  the.  data  and  the  state.  The 
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development  here  will  consider  first  the  restricted  case  where  the 
data  only  is  affected  (e.g. ,  a  random  component  of  station  location 
error)  since  most  of  the  important  principles  are  illustrated  with 
this  simpler  model.  The  technique  will  then  be  extended  to  the  more 
general  case.  We  begin  by  stating  the  familiar  data  equation 

*  Ax  +  By^  +  n^,  (5.11) 

where  x  is  the  state  being  estimated  (dimension  n) 

y^  is  the  instantuncous  value  of  the  consider  parameters 
(dimension  p)  and  the  data  has  been  pre- weighted  so  that  the  additive 
noise  n^  has  zero  mean  and  unity  covariance.  The  transition  equation 
for  the  y's  is  given  by 

yi+l  E  +  vi  (5-12) 

where  |  |M|  |<  1  to  keep  the  process  from  growing  without  bound.  *  The 

w's  have  the  same  dimension  as  y^,  and  are  generated  as  white  sequence 
noise  having  zero  mean.  Normally  M  and  the  covariance  on  Wj  would  be 
chosen  so  that  the  covariance  of  y^  would  remain  stationary. 

In  terms  of  these  quantities,  the  problem  is  v  valuate  the 
performance  of  the  filter  which  estimates  x  but  ignores  the  influence 
of  thf  y  variables.  To  this  end,  suppose  that  at  time  t^  just  prior 
to  the  data  taking  at  that  point,  nn  information  array  is  obtained 
that  contains  the  current  estimate  of  x  and  the  true  and  computed  co- 
variance  of  x,  y^,  and  r,  where  p;  is  an  n  dimensional,  zero  meen, 
random  variable  of  unity  covariance  whose  involvement  in  the  problem 
will  become  clear  subsequently. 
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(5.13) 


From  (5.13),  the  estimate  is  given  by 

x  -  Rx_1  z[  (5-14) 

The  computed  (i.e. ,  the  stochastic  process  free)  covariance  is 

Ai  ■  V1  V11  <5.15) 

c 


and  the  true,  or  consider  covariance  is: 


(5.16) 


Equation  (5*16)  can  be  viewed  in  either  of  two  ways:  l)  It  arises 
from  the  concept  of  the  consider  option  and  is  simply  an  application 
of  (5.3)  by  recognizing  that  the 


is  the  inverse  square-root  on  the  covariance  of  the  consider  parameters. 
2)  Because  of  the  structure  of  the  right-hand-side  of  5.13>  even  if 
the  consider  parameters  were  estimated,  their  estimates  would  be  zero. 
Thus,  the  left-hand-side  of  (5*13)  can  be  considered  as  the  inverse 
square  root  of  the  covariance  of  all  the  parameters.  In  that  sense 
(5.16)  represents  the  upper  left  nxn  portion  of  the  full  covariance. 

A* 

7* 


Data  Processing 


Data  is  taken  at  t^  in  accordance  with  (5.11),  and  added  to  the 
existing  dnta  in  the  following  manner: 
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Equation(5>17)  is  compacted  with  a  standard  Householder  transformation 
resulting  in  the  first  n  rows  as 


[R*  ,  R'  ,  R» 
x  *  x^  *  xy 


A 

X 


(5.18) 


This  data  insertion  operation  has  been  performed  as  if  nothing  below 
the  first  row  existed.  This  is  best  justified  on  a  simpler  problem. 
Suppose  there  are  but  two  classes  of  variables,  x  and  y  as  in  (5.1). 
That  is. 
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T 

X 

xy 
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0 

Ry  _ 

y 

0 

-1  -IT 

wheze  R  R  is  the  covariance  on  y.  Additional  data  is  taken  and 

y  y 

adjoined  as  follows. 


Row  interchanges  can  be  made  without  changing  the  problem.  Re-write 
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the  above  as 


A  Householder  transformation  is  applied,  resulting  in 


Notice  that  R'  ,  R'  obtained  in  this  fashion  is  identical  to  what 
x  xy 

would  have  been  obtained  if  the  transformation  had  been  applied  to 


Further,  R'^  is  not  needed.  The  covariance  on  x  is  now 

(R’"1)(R'  "1)T  +  R'"1  R'  A  R'  T  (R'  ~1)T 
vx''x  x  xy  y  xy  x 

This  completes  the  justification. 


Returning  to  the  full  problem,  the  estimate  at  this  juncture  would  be 
x  -  R  x  *1 

and  the  computed  and  consider  covariance  would  be  calculated  as  in 
(5.15)  and  (£.16).  The  primed  quantities  would  be  used,  of  course; 


the  covariance  on  c  and  y^  would  remain  unchanged. 
Mapping 


Equation  (5*18)  must  be  referred  to  t^+1>  This  is  done  by  solv¬ 
ing  (5.12)  for  yit 
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yi  "  M’lyi+i  "  M"lvi 


and  re-writing  5.18  in  terms  of  x,  yi+^»  end  We  obtain 


"  R'  R'  R' 

x  x$  .  xy 

0  I  0  . 

E 

-1  .  -1 

0  0  Ry 

M  yJ+1  -M  wt 

(5.19) 


Performing  the  indicated  matrix  multiplication,  effecting  certain  row 
interchanges  and  adding  the  a  priori  information  on  w: 


(5.20) 


R'  R’  -R'  M"1  R*  m"1 

x  x^  xy  xy 

k 

0  10  0 

5 

0  0  -R  M-1  R  M-1 

y  y 

wi 

0  0^0 

■ - 1 

>r 

_ i 

The  last  row  of  5.20  is  to  reflect  that  w^  has  mean  zero  and  covari¬ 
ance  /^.  The  most  convenient  way  of  expressing  this  is  to  add  a  p 
dimensional  dota  vector  whose  a  priori  value  is  zero  related  to  w 
through  an  identity  transformation.*  That  ic 

Zs  -  O'  -  IWj  +  «„ 


ET.J  .  0  and  e[V„T] 


The  dummy  measurements  are  weighted  by  pre-multiplying  by  to 

obtain  i 

0  ■  K  'li*% 


where  Efn  n.  T1  -  I 
w  w  j 


♦This  particular  operation  was  first  used  by  Dyer  and  McReynolds  in 
their  derivation  of  an  optimal  square-root  sequential  filter  for 
process  noise  (Reference  12). 
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which  ia  the  lust  row  ol'  (5.20). 

If  (5*20)  is  compacted  with  a  Householder  transformation,  the 
following  array  is  obtained. 


where  the  's  have  been  retained  to  emphasize  that  the  1st  2  rows  re¬ 
main  unchanged  by  this  trans format ion.  The  lower  right  2x2  portion 
represents  the  inverse  square-root  of  the  covariance  ow  w  and  y^+1* 
These  are  correlated  stochastic  variables.  If  the  contribution  of 


their  effects  to  the  true  covariance  were  calculated  at  this  juncture, 
it  would  be 


For  reasons  to  be  clear  presently,  it  is  desired  to  replace  the 
^i+1*  Pa*r  a  new  Pft*r>  ^i+l*  w^ere  f  is  a  p  dimensional  random 
variable  independent  of  y^+^  (and  5)  with  covariance  (for  convenience) 
I.  In  other  words,  abandoning  the  matrix  notation  for  the  moment  and 
writing  out  the  first  row  of  (5.21) 


7,,.  »  R*  x  +  R  *  +  R  w.  +  R  y.  +  n 

i+1  x  x$  xw  i  xy  *i+l 

We  wish  an  equivalent  data  vector  of  the  form 

zi+l  "  R'xx+RV  +R^?  +  R,xyyi+l  +  n 


(5.22) 

(5.23) 
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For  equivalency,  two  conditions  must  be  observed: 


1)  the  contribution  of  ~  and  y^+^ 

to  the  total  error  covariance  from  5*23  must  be  identical  to  the  con¬ 
tribution  of  vi  and  y^+^  from  5.22). 

2)  the  eras:;  correlation  of  with  y^+^ 

must  be  the  seme  in  5*23  and  5.22.  This  is  needed  when  adding  new 
data  which  involves  yi+^.  Since  in  effect  5*23  ploys  the  role  now  of 
apriori  data,  its  relationship  to  y^+^  must  be  expressed  properly. 
Invoking  the  2nd  condition  first, 

*n,  E  [VmT]  +  V  E  [ymJ,l+lT]  -  »'xy  E  W 
Recoiling  that 

yiU  ■  ,<yi  +  ”l 
•  then  Er«i  yltlT]  " 


end  h  yi+1  y1+1 


U  h-1b-1t 


(in  the  usual  stationary  circumstances,  Ay  will  not  be  a  function  of  i) 


Thus, 


R'  b  R  A  A  R 

xy  xw  w  ‘  y  xy 

the  first  requirement  implies  that 
lKxw  RxyJ  F  Aw  Awl  T  ^ 


(5.24) 


l> 

«  Fr  ,  r»  '!  f  i 
l.  h 


(5.25) 


°  ir*Txr 


Carrying  out  the  multiplication  indicated  by  5.25»  substituting  for 
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value  obtained  for  R'xy  in  (5.2U),  and  solving  for  Rx~ 


T 

R  ~  R  ■*  R  A 
x$  Xf  xw  'H i 


xy 

T  -IT 

R  ‘  K  -  N,  Ay  Aw  R 


xw 


xw 


XV/ 


(5.26) 


R*,  Aw  (1  -  V'1  V'T)  \J 


Finally  since  a  5. 25  can  be  written  in  explicit  square  root 
form  as 

£  ft  .  *  7  a  ■A  *  (5.27) 


Rx?  -  R*v  (I  -  ij  V 1  AwiT)i 


Since  r  and  “g  are  independent  of  any  of  the  other  variables ,  their 
influence  can  be  broken  out  separately  (this  of  course  is  the  motiva¬ 
tion  foi  the  transformation  to  r,  y^+^  fro®  yi+l^*  This  contribu¬ 


tion  is  written  as 


R  “1  I  R  ,  R  ~J 

X  L  XJ*  XJJ 


xc 

RT~ 


-IT 


(5.28) 


Now,  if  P  P  «  I,  an  orthonormal  transformation  is  inserted  in  the 
middle  of  (5.28)  and  PT  is  such  that 


XS 


XS 


H 


Then,  the  contribution  from  g  and  g  can  be  written  in  terms  of  a 
single  random  variable  again,  which  shall  continue  to  be  denoted  as  *. 


Finallyjthe  post  mapping  information  array  is  written  as 

R  •  K  R'  n  r 

x  x*  xy 


0 

0 


I 

0 


y  J 


x 
§ 

yi+u 


'i+l 


0 

0 
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But  equation  5.29  is  of  the  same  tom  as  5.13  except  that  one  full  cycle 
of  data  taking  and  mapping  has  been  completed.  The  reason  for  the 
original  introduction  of  the  5  variable  should  now  be  clear.  It  is  to 
express  the  errors  in  x  that  arise  from  accumulated  effects  which  are 
orthogonal  to  the  current  y.  Its  dimension  is  the  some  as  x. 

Extension  to  the  General  Case 

If  it  is  now  desired  to  consider  the  general  case  where  the  pro¬ 
cess  affects  not  only  the  data  but  the  state  as  well,  the  procedure 
differs  considerably  in  detail  but  not  in  concept  from  that  just  des¬ 
cribed.  When  employing  a  sequential  estimation  that  optimally  treats 
process  noise,  no  augmentation  of  state  is  necessary  since  the  indirect 
effect,  of  the  process  noise  on  the  data  through  the  perturbation  of 
state  is  automatically  accounted  for  when  the  covariance  matrix  of 
state  estimation  errors  is  increased  by  the  process  noise  during  the 
mapping  operation.  In  the  case  where  we  are  considering  the  process 
noise,  the  information  array,  R^,  expresses  the  information  concerning 
a  hypothetical  process  noise  free  system.  An  auxiliary  variable  is 
needed  to  play  the  role  of  representing  the  accumulated  effects  of  the 
process  noise  on  the  current  data.  Define 

x(0)  »  initial  state  of  the  system  at  tQ 
Xj,  ^  current  system  state  at  tj 

^x^(O)  £  current  effective  departure  from  the  original  initial 
system  state  due  to  process  noise.  ax^(°)  is  such  that 

*i  -  U  (tif  tQ)  (x(o)  +  Axt)  (5.29) 

Figure  5.1  illustrates  the  meaning  of  Axj  pictorial ly.  Thus,  for 
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STATE, 


Figure  5.1.  Definition  of  Ax  (t±) . 
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convenience  we  shall  assume  we  are  estimating  not  the  current  state, 
but  the  original  state  x(o)  and  are  considering  ^  find  y^.  From 
Figure  5.1,  if 

xi+l  ’  u  (t1+l»  \)  +  Gi  yL  ( 5# 30) 

then 

Axni  B  Axi  +  u  ^i+i*  to)  Giyi 


Considered  alternatively,  the  ^x^,  y^  jointly  now  play  the  role  as  did 
y^  alone  in  the  previous  development.  Their  combined  state  transition 
would  be 


—  — 

..-1  " 

—  - 

Ax 

B 

I 

u  0 

Ax 

+ 

0 

y 

0 

M 

y 

1 

L-  J 

i+1 

- 

- 

i-  -J 

i 

-  - 

(5.31) 


The  data  taken  at  t^  will  relate  to  the  instantaneous  state  x^  (or 
through  5.29  the  sum  ol’  x(0)and  end  the  current  process  noise, y^. 


«  A(x(0)  +  )  +  B  y^  + 

Zi 

■  A  x(0)  +  A&x^  +  B  y^  +  n^ 


(5.32) 


The  information  array  just  prior  to  data  taking  at  t^  (analogous  to 

5.13)  is 


(5.33) 


R 

R 

R 

R 

x(0)  " 

~  zi" 

X 

x* 

x/vx 

xy 

0 

I 

0 

0 

% 

m 

0 

0 

R 

R 

/Xj 

tx 

/xy 

0 

0 

0 

0 

Ry 

yi 

- 

- 

- 

- 

t  A 
L 


0 


Z  I 

!  i  1 
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where  the  last  row  shows  the  subsequent  addition  of  the  data  at  t^ 
from  5»32. 

The  procedure  follows  that  used  previously  with  estimation  and  co- 
variance  computations  analogous  to  5*l4  through  5*16.  After  data  is 
added  at  t^,  compaction  takes  place  as  before  (see  5*17  and  5.18).  In 


order  to  map,  5*31  is  solved  for 


obtaining 


-i  -1~ 

“  ■*“ 

-i  -i 

Ax 

m 

I 

-U  GM  A 

AX 

+ 

U  GM 

y 

0 

M-1 

y 

■m”1 

*-  — i 

i 

- 

- 

- 

i+i 

L  J 

(5.34) 


This  is  substituted  in  5*33  after  compaction,  the  same  algebraic  manf 
pulation  takes  place,  the  a  priori  data  on  w^  is  inserted,  and,  once 
again,  a  forward  triangularization  is  performed.  The  following  array 
is  obtained: 


R  R  R  R  R 

x  xf  xw  x^x  xy 

x(0) 

0  10  0  0 

? 

0  0  R  R  .  R 

w  wax  wy 

Wi 

0  0  0  R  R 

AX  A*y 

A*i+1 

0  0  0  0  R 

y 

yi+l 

(5.35) 


Again  an  .independent  rendorn  veriable  •  must  be  introduced  to  replace 

IV} 

w  .  The  Ploy  the  same  role  as  did  y.  .  before.  Still,  the 

i  i~y  J  i+1  i4l 

algebra  is  enough  different  to  present  that  transformation  ir.  some 


detail. 


““['L 


L 

B 


0 
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(5.36) 


and 

We  wish  to  write 


r  R  ,  P  ]  -  R 
-  x^x*  xyJ  xq 


z  ■  R  x  +  R  p+R  w  +.R  _  o 

X  X£  ■  XW  Xft  w 


as  an  equivalent 
z 


K  X  +  Kw  ?  +  Ryr  W  +  R'xn  0 
x  x?  *  xo 


For  equivalency  it  is  required  that 

"xw  n1]  +  R*o  E[°  °T]  ’  R'*n  E[n  nT] 

Defining 


*[>0*]4  avo 
*[n  nT]  *. ^ 


We  obtain 


R'  -  R 
Xfi  xw 


WAn 


-l 


+  R 


xo 


We  now  must  select  R  ~  r.o  that  p  and  ft  will  contribute  ar. 
variance  in  (5*37)  as  do  w  and  ft  in  (5.36).  This  is  done 
(5.26)  through  (5.27).  In  this  case  the  result  is 


R  rw  R  n  *  R  A 
Xf  X?  xw  ^ 

It  is  easily  shown  that 
**»  ■  [°  *»] 


-1 


xw 


Rxw  ^Wft  ■  ^W ft  Rxw 


and 


*0 


-I 


R 


,  X 

R 

R  1 

R 

Ax 

Ax 

Ax 

■t*  y 

T 

It 

R 

r  1 

R  + 

Axy 

A* 

Axy 

Axy 

T 

R  A  R 
y  -3 


(5.37) 


(5.38) 


(5.39) 

much  to  the 
just  as  in 

(5.*»0) 
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Finally 


R  rv  "  R 

Xf  XV 


and  aa  before, 
procedure . 


*•*  {  1  -  (*.**>  K*y  Vy  +  Rly  Ry]  *.* }  *  (5 M) 

R  and  R  ~  are  combined  completing  the  generalized 


Propagating  Forward  and  Evaluation 

The  square-root  of  the  information  matrix  (5 - 33)  can  be  referred 
to  a  new  epoch  at  any  time  by  simply  post  multiplying  by  the  inverse 
state  transition  mairix,  i.e., 


‘  R  U_1 

R 

R  .  if1 

R 

x(j)  " 

z 

x 

XS 

xAx 

xy 

0 

l 

0 

0 

5 

0 

0 

0 

V"'1 

Vp 

a^U) 

m 

0 

0 

0 

0 

R 

p 

0 

P 

_+ 

where 

U  «  U  (t  y  to) 

and 

(d)  ■  U  (*j,  tQ)  AXj 

In  evaluating  for  the  covariance  of  the  error  in  the  estimate,  it 
must  be  remembered  that  the  error  in  the  estimate  of  x(0)  is  the  error 
of  the  original  injection  conditions.  It  is  usually  desired  to  know 
the  error  in  the  effective  conditions,  i.e., 
x  -  x  ♦  (ax  -  L$t) 

A 

Since  Ax  is  always  zero,  we  equivalently  wish  to  calculate  the 
E  T  { (x-x)  ♦  Ax|{(x-£)  +  Ax  l  ] 
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It  is  easily  shown  that  this  is 


X£X 


R  “1T  +  R 
x  x 


-1 


X£X 


A  R  T  R 


-IT 


where 


>x 


R  "1  (R  “1)T 
AX  fix 


and  is  the  covariance  of  perturbations  due  to  the  process  noise. 

It  may  appear  somewhat  unorthodox  to  have  cast  this  problem  in 
terms  of  estimating  the  initial  state,  particularly  in  a  problem  that 
seems  most  naturally  to  lend  itself  to  estimating  the  current  state 
in  the  usual  sequential  manner.  This  has  teen  done  for  two  reasons: 


1.  The  algorithms  so  derived  will  mate  more  easily  with 
already  existing  computational  machinery  for  their 
application  in  Chapter  VI. 

2.  It  should  be  clear  that  with  minor  modifications  (all 
simplifications)  to  whot  was  presented,  we  could  have 
produced  an  algorithm  for  optimally  estimating  x(0) 
and  ax< .  That  is,  a  square-root,  one-pass,  forward 
smoother  in  discrete  form  would  have  resulted. 

Forward  smoothers  are  discussed  by  Nishimura,  Ref. l6. 


A  Numerical  Example 

To  clarify  the  procedure  and  to  demonstrate  its  validity,  a 
simple  numerical  example  was  constructed,  one  where  the  proper  answers 
could  be  calculated  by  appeal  to  a  much  simpler  point  of  view.  Suppose 
both  x,  w  and  hence  ax  have  dimension  one.  Further,  let  y  *  w,  i.e, 
there  really  is  no  need  to  deal  with  the  auxiliary  variable,  y. 

Finally  let  the  measurements  directly  measure  x  and  ax: 

Zj  ■  x  +  ax^  +  n^ 


where 


AXt 


*Vl  +  wi-l 


and 

M  ’  °*  EInt  "j]  •  1  1  ■  3 

0  1  *  6 

ErAXi J  «  0,  if  AXl2]  -  1 

E[*i]  .  0,  e[„4  wj]  1  -  i 

o  i  i  j 

A  sequence  of  measurements  is  obtained: 
zl  -  x  +  gx^  ♦ 

*2  -  x  +  &2  +  "2  “  x  +  A*i  *  wi  4  n2 
*2  “  x  ^  ^2  +  **3 

L 


\  ’  *  +  *X1  *  *  V1  +  ”n 
i»l 


The  filter  is  aware  only  of  the  presence  of  x  and  n.  Thus  it  will 
have  the  form: 

x  »  jj  S  *t 

"  i-1  1 

The  true  performance  can  be  calculated  by  considering 


A*!  4  r  wt  ♦  n 
x  1-1  1  n 

to  be  the  true  measurement  noise  and  calculate  the  actual  filter  per- 
fomance  which  equally  wights  tbs  measurement*.  It  is  easy  to  show 
that  for  aa  equal  weighting  filter  where  the  actual  data  noise  it 
tha  true  covariance  of  estimation  srror  is 


displays  results  from  crch  method  for  the  first  fev  data  points.  It 
is  seen  that  acrecncnt  is  maintained  to  8  place*  accuracy,  the  printout 
resolution  used  (computation  vas  performed  in  double  precision). 
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Table  5.1 


Stochastic  Consider  Option  Check 

No.  of  Data  Points  Analytical  Result  Stochastic  Consider 

________________  _  Option 

Result 


1 

2 

2.0000000 

2 

7/4  -  1.75 

1.7500000 

3 

17/9  -  1.6888888. . 

I.8888889 

4 

34/l6  *  2.125 

2.1250000 

5 

60/25-  2.4 

2.4000000 

as 


CHAPTER  VI. 


THE  LONG-TERM  INFLUENCE  OF  POORLY  MODELED  FORCES  AND 
STATION  LOCATION  MOVEMENT  ON  ESTIMATION  ACCURACY 

In  Chapter  IV,  the  short-term  analytic  doppler  and  range  track¬ 
ing  models  were  applied  to  the  problem  of  determining  the  influence 
of  tracking  station  location  errors  end  non-gravitational  forces  on 
state  estimation  accuracy.  In  particular,  the  non-gravitational 
force  problem  was  explored  quantitatively  for  1)  several  trajectories^ 
2)  Inclusion  of  range  data  or  not,  and  3)  attempted  solution  for  the 
unknown  force  or  not.  The  thread  of  that  investigation  is  picked  up 
here  and  both  the  station  location  problea  and  the  small  force  problem 
are  explored  quantitat  vely,  this  time  with  more  precise  modeling  and 
over  long  tracking  intervals. 

In  this  Chapter,  vc  shall  try  to  take  an  operational  point  of 
view  and  address  problems  inherent  in  using  tracking  data  during  a 
real-time  mission  when  the  objectives  arc  to  best  determine  the  state 
and  predict  the  encounter  conditions  of  the  trajectory  relative  to 
the  target  planet.  This  la  a  problem  that  is  distinct  from  the  post- 
flight  data  reduction  problem  which  emphasises  careful  solutions  of 
astrodynamlcai  constants,  tracking  station  part:  •.*  *  calibrations,  snd 
perhaps  refined  trajectory  estimates. 

Missions  up  through  Mariner  *Cy  have  relied  upon  very  simple 
filters  in  the  sense  thnt  the  parameter  list  if  abbreviated,  usually 
only  includliv!  the  six  components  of  trajectory  state  and  solar 


®9 


pressure.  Tracking  data  from  past  missions  is  used  to  produce  the 
most  careful  calibration  possible  as  to  the  station  locations  and  the 
important  astrodynomical  constants  that  affect  the  solution  using  the 
abbreviated  filter.  During  the  flight,  Universal  Tima  and  Polar  Motfcn 
which  are  essential  to  locate  the  instantaneous  position  of  the  sta- 

m 

tions  relative  to  an  assumed  known  inertial  coordinate  frame,  are  mon¬ 
itored  closely.  Depending  on  the  accuracy  required,  the  troposphere, 
ionosphere,  and  space  plasma  -  oil  factors  directly  affecting  the  data 
•are  calibrated  in  real-time  by  means  external  to  the  orbit  determina¬ 
tion  process.  Finally  conventional  wisdom  has  it  that,  for  example 
in  predicting  encounter  conditions,  only  a  relatively  short  arc  of 
data  should  be  used  in  the  final  estimation  process.  This  data  arc 
should  not  exceed  30  days  and  would  normally  be  the  last  data  received. 

The  reasons  for  the  abbreviated  solution  vector  and  the  short 
tracking  arc  used  can  be  illuminated  here  but  briefly  (for  a  more 
detailed  discussion  of  these  factors  and  principles  as  applied  to 
planetary  cruise  orbit  determination,  see  Reference  17): 

1.  Station  location  errors,  even  carefully  calibrated  as 
they  are,  still  are  the  predominant  error  source.* 

Further  it  is  easy  to  demonstrate  that  by  using  ?  or 
3  conths  of  data  prior  to  encounter  and  estimating 
the  station  coordinates  produces  a  claioed  orbit  de¬ 
termination  performance  substantially  superior  to 
that  of  using  the  shorter  tracking  arc  and  consider¬ 
ing  the  station  errors  at  their  npriori  valuer.  Vhy 
not  then  follow  this  apparently  more  optimal  procedure? 

The  first  flaw  is  that  thin  same  long  arc  analysis  will 
also  reveal  that  the  station  longitudes  will  be  deter¬ 
mined  to  better  than  .'j  a,  standard  deviation.  This 
is  not  reasonable.  We  have  been  using  nearly  a  decade 
of  Interplanetary  tracking  data  and  heve  never  obtained 


•Currently  the  DSil  trreking  station  locations  are  known  approximately 
to  1.5  m («  )  in  distance  off  the  spin  axis  and  3  c(~  )  in  the  lonM* 
tude  dl-  rs  reetSc.n  r> 
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consistency  of  answers  much  better  than  the  claimed  know¬ 
ledge  of  3  m.  It  is  not  likely  that  an  additional  months 
or  so  of  tracking  will  markedly  improve  this  situation. 
Therefore,  the  long  arc  analysis  results  are  optimistic, 
and  it  is  not  clear  that  the  actual  result  would  in  fact 
not  be  worse  than  the  simpler  short  arc  result. 

2.  The  small  unmodelcd  forces  may  be  producing  actual  filter 
divergence  when  too  long  an  arc  of  data  is  used.  Filter 
divergence  is  often  characteristic  of  mismodeling  of  dynamic 
parameters,  but  until  recently,  there  has  been  little  hard 
evidence  as  to  whether  this  is  the  case  for  interplanetary 
orbit  determinations  over  the  time  span  we  are  interested 
in.  In  the  face  of  thi3  paucity  of  evidence,  divergence 
has  remained  a  vague  fear  and  prudence  has  dictated  that 
short  tracking  arcs  be  used.  Chapter  IV  has  concluded  that, 
because  the  data  employed  arc  almost  directly  affected  by 
accelerations,  unmodeled  forces  immediately  produce  errone¬ 
ous  results  at  the  levels  determined,  but  this  analysis 
gives  no  clue  as  the  their  effects  in  the  longer  term. 
Recently  J.  F.  Jordan  et  al,  produced  results  which  indi¬ 
cated  that  divergence  did  not  occur  over  several  months 
tracking  intervals  when  constant  forces  were  simply  ignored 
in  the  filter  itself  (Reference  18). 

We  hope  to  shed  some  light  on  both  station  location  and  unmodeled 
force  effects  here  as  an  application  of  the  computational  machinery 
developed  in  Chapter  V. 

Station  Location  Error  Study 

The  reason,  of  course,  that  the  formal  statistics  for  a  long 
data  arc  including  station  location  solutions  are  optimistic  is  be¬ 
cause  of  too  simplistic  a  model.  In  solving  for  station  parameters, 
the  difference  between  the  actual  station  location  and  that  used  in 
the  abstract  model  is  assumed  to  be  a  constant  vector  in  the  Earth's 
rotating  co-ordinate  frame.  In  actuality,  this  difference,  and  hence 
the  station  location  error,  has  a  significant  stochastic  component. 
This  component  arises  from  temporal  variations 

i.  Universal  Time  calibration  error 
ii.  Polar  Uotlor.  cr!1br*.»! -n  '*rr:»r. 


in  the 
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In  addition,  uncalibrated  tidal  effects  may  be  producing  signi¬ 
ficant  station  location  errors.  Finally,  the  uncalibratcd  portion  of 
the  ionosphere  and  space  plasma  tend  to  produce  unmodeled  diurnal 
signatures  in  the  data  that  arc  effectively  station  location  shifts. 

There  is  little  hope  that  these  effects  can  be  modeled  as  Marko¬ 
vian  stochastic  processes  accurately  enough  to  include  them  as  a  form¬ 
al  part  of  the  filter  model.  It  should  be  possible  to  postulate  rea¬ 
sonable  enough  models  for  these  processes  to  be  used  as  consider 
variables  however.  For  example,  suppose  we  wish  to  test  two  rather 
simple  filters: 

1.  Estimate  state  only 

2.  Estimate  state  and  a  constant  station  location  error. 

These  two  filters  can  be  evaluated  in  the  presence  of  both'  the 

fixed  and  stochastic  station  errors  for  a  range  of  possible  descrip¬ 
tions  that  is  thought  reasonable  for  the  stochastic  portion.  If 
filter  2  outperforms  filter  1  over  this  whole  range,  confidence  will 
be  gained  that  it  truly  is  superior.  Moreover,  at  the  same  time 
statistics  would  be  generated  for  filter  2  that  could  be  believed. 

In  short,  the  Bax  lx  operating  here  is:  The  filter  evaluation 
should  include  a  model  that  is  at  least  one  step  (in  some  sense)  more 
sophisticated  than  the  filter  Itself. 

In  using  filter  1),  that  step  has  been  the  consider  option  in¬ 
cluding  the  static  station  location  errors.  The  results  can  be  be¬ 
lieved  and  probably  would  not  be  much  worse  (they  night  be  better)  if 
the  refined  stochastic  point  of  view  were  utilised.  In  using  filter 
2),  a  stochastic  portion  rust  be  postulated  if  the  statistical 
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results  are  not  to  be  hopelessly  optimistic. 

Before  proceeding  to  discuss  the  results ,  it  should  be  emphasized 
that  while  every  attempt  has  been  made  to  use  realistic  numbers,  a 
formal  and  authoritative  model  of  the  stochastic  station  location  com¬ 
ponents  is  beyond  the  scope  of  this  work.  We  are  pursuing  a  point  of 
view  on  technique  here  and  the  results  can  at  best  be  claimed  as  rep¬ 
resentative.  Performance  of  the  actual  studies  which  bring  change  to 
current  operational  philosophy  are  not  best  carried  out  in  the  rela¬ 
tive  isolation  present  when  pursuing  independent  research. 

To  obtain  the  desired  numerical  results,  a  Univac  1108  computer 
program  was  designed  and  built  that  executed  the  Generalized  Consider 
algorithms  developed  in  Chapter  V.  Since  it  was  desired  to  concen¬ 
trate  on  the  estimation  aspects  of  the  problem,  no  trajectory  or  par¬ 
tial  derivative  packages  were  built.  Rather,  the  program  was  designed 
in  such  a  way  that  it  could  be  driven  by  JFL's  Double  Precision  Orbit 
Determination  Program  (Ref.  19)  and  simply  use  the  trajectory  informa¬ 
tion,  state  transition  matrices,  and  data  partials  generated  there. 

In  this  way,  we  become  heir  to  all  the  sophisticated  modeling  and 
generality  built  into  tho  DFODP  with  no  effort  on  our  part. 

The  point  of  view  taken  on  the  stochastic  portion  of  the  station 
locations  (and  in  the  following  sub-section  •  on  the  random  forces) 
was  that  they  were  piecewise  constant  over  an  Interval  £t.  Then  from 
Interval  to  interval  they  were  correlated  as 

#-dt/T 

As  long  as  At  «  ▼,  this  is  thought  a  reasonable  approximation  of  s 
continuous  exponentially  correlated  s tocher  tie  process.  With  thtr 


approximation,  the  data  partials  could  be  processed  in  batches  with 
At  time  span  and  greatly  reduce  the  computations  needed  compared  to 
that  required  if  all  of  the  steps  of  the  stochastic  consider  option 
were  executed  each  data  point.  All  of  the  station  location  results 
presented  here  assume  a  correlation  time,  t,  of  25  days  and  use  a  At 
interval  of  U  days.  Test  results  using  lower  At’s  yielded  results 
only  negligibly  different  from  those  shown  here. 

The  Specific  Case 

Again  a  Viking  *75  Mars  trajectory  was  selected  but  this  time, 
the  simulated  data  commenced  six  months  prior  to  encounter  so  that 
arbitrarily  *>T,/.s  of  data  could  be  investigated.  One  pass  of 
doppler  only  data  from  a  single  station  each  day  was  used  at  an  effec¬ 
tive  data  noise  of  1  mm/sec  for  a  1  minute  sample.  The  station  was 
assumed  to  have  station  location  errors  as  follows: 


Constant 

Stochastic 


Figures  6. 1-6. 3  display  the  performance  of  the  filter  which  cs- 
tlaates  the  6  components  of  state  only.  Shown  are  the  uncertainties 
of  the  three  components  of  initial  position  in  the  plane-of-sky  co¬ 
ordinate  system  used  in  Chapters  III  and  IV: 
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Figure  6,1  Ror.ce  Coeponcnt  Filter  Performance  State  Only  Versus 
Tracking  Tine  -  Doppler  Only. 
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Figure  6.2.  Declination  Component  Filter  Farfomance  State  Only 
venue  Tracking  Tine  •  Doppler  Only. 
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Figure  6.3.  Right  Ascension  Conponents  Filter  Performance  State 
Only  versus  Tracking  Time  -  Doppler  Only. 
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Earth-probe  range  -  r 
Geocentric  Declination  -  ft 
Geocentric  Right  Ascension  -  a 
There  arc  three  curves  shown  in  each  figure: 

1.  Standard  deviation  that  the  filter  itself  computes  (lowest 
cuive) 

2.  Standard  deviation  assuming  only  the  existence  of  the  con¬ 
stant  station  location  errors  (middle  curve) 

3*  Standard  deviation  assuming  both  the  existence  of  the  con¬ 
stant  end  the  stochastic  portions  of  the  station  location 
errors.  This  curve  represents  the  true  filter  performance 
for  the  complete  physical  mode)  assumed  (highest  curve) 

It  can  be  seen  that  the  second  curve  is  a  reasonable  representa¬ 
tion  of  the  true  performance.  That  is,  if  the  filter  ignores  the  sta¬ 
tion  location  errors,  it  is  not  critical  to  recognize  that  they  have  a 
stochastic  component  when  trying  to  evaluate  that  filter's  performance. 
The  lowest  curve,  of  course,  is  not  a  good  measure  of  performance. 

Figures  6.U  and  6.5  plot  similar  information  for  the  filter  that 
estimates  the  constant  station  location  errors  but  is  unaware  of  the 

existence  of  the  stochastic  component.  Again,  the  lower  curve  in  each 

■ 

figure  represents  the  filter's  own  computation  of  its  performance 
which,  this  time,  includes  the  constant  station  location  error  model. 
Nov  it  is  critical  to  include  the  stochastic  portion  in  the  evaluation 
in  order  to  get  reasonable  results.  The  true  performance  is  given  by 

r  i  I  . 

the  top  curve. 

Cross-comparison  of  the  top  curves  for  the  two  filters  indicates 
that  ei^imating  the  station  location  errors  does  refine  the  results; 
it  does  not  however,  revolutionize  them.  Table  6.1  facilitates  this 
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Figure  6.4.  Declination  Component  Filter  Performance  State  and 
Station  Location  Solution  vereuc  Tracking  Time  - 
Doppler  Only. 
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TRACKING  TIME,  doyi 

Figure  6.5.  Right  Ascension  Component  Filter  Performance  State  and 
Station  Location  Solution  versus  Tracking  Time  - 
Dopf  ler  Only. 
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cross-comparison  by  listing  the  actual  standard  deviations  for  all 
sin  of  tha  co-ordinates  after  126  days  of  tracking  data.  Although 
•ach  co-ordinate  haa  been  Improved  by  extending  the  parameter  act, 
it  la  nowhere  near  aa  draoat'c  aa  if  the  atochaatlc  portion  were  not 
preaent.  Significantly  absent  la  a  large  improvcawnt  in  the  declina¬ 
tion  direction,  which  la  always  a  major  contributor  to  overall  per¬ 
formance  errors.  The  per  cent  improvement  in  a  and  the  three  veloci¬ 
ties  la  more  significant  -  better  than  a  factor  of  two  in  each  case. 

Finally,  ve  wish  to  preaent  some  data  coooe ruing  the  accuracy 
and  validity  region  of  the  analytic  model  discussed  in  Chapter  IV. 

This  is  found  in  Table  21  which  converts  the  effect  of  the  eons  tent 
station  location  errors  computed  analytically  with  those  obtained  from 
the  BfODP  vhoee  model  aaaumaa  little  else  than  teuton  and  Klasteln 
were  correct.  Comparing  the  "total  no"  line  with  tha  8  day  result, 
remarkable  agreement  is  found  in  all  but  the  r  position  sad  declina¬ 
tion  velocity  components.  Kven  hers,  the  agreement  is  quite  respec¬ 
table.  As  more  tracking  data  is  taken  (and  assertions  in  the  snaly- 

.  4 

tlo  model  are  violated),  the  agreement  deteriorates  aa  would  be 
expected.  Of  particular  note  is  that  right  ascension  position  and 
velocity  estimates  become  less  sensitive  to  station  location  errors 
as  more  data  is  included.  The  increase  in  the  range  error  with  tins 
could  probably  be  removed  with  the  addition  of  ranging  data,  although 
this  must  be  carefully  done  lest  an  even  greater  error  appear  in  some 
other  co-ordinate. 

Small-Force  Study 

A  study  of  non-nravitatlonal  forces  such  aa  presented  in 
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Table  *.2 

Corperlson  of  Chapter  XV  Analytic  Model  with  Exact  Model 


b02xl0*7  h.26alO*7  .19*10' 


Chapter  IV  will  provide  at  rone  Motivation  to  try  to  incorporate  then 
somehow  in  tha  estimation  a c hear.  Juat  how  to  do  this  it  not  clear, 
however.  The  Main  difficulty  la  that  not  only  art  tha  forcaa  aa  a 
function  of  time  unknown,  vary  little  information  la  uaually  available 
to  even  describe  the  a  toe  has  tic  proceaa  from  which  they  are  drawn. 

Tha  auln  sources  of  tha  force  a  are: 


1.  lho  aolar  pressure  force  and  lta  interaction  with  the  space- 
craft.  To  a  first  approximation  this  can  be  aodeled  aa  a 

1 Jr  force,  nuch  Ilka  aolar  gravity  itself,  except  it  is  non- 
conservative.  Depend  inc  on  the  axial  symmetry  of  tha  spac¬ 
ers!*,  this  force  will  be  principally  directed  outward  from 
the  sun.  In  practice  tha  tain  source  or  uncertainty  here  is 
in  predicting  tha  affective  albedo  of  the  spacecraft  (parti¬ 
cles  that  are  reflected  will  1  apart  twice  the  momentum  aa 
those  absorbed)  and  lta  time  varying  characteristics  (  as 
the  spocecralt  ages,  the  reflective  properties  change,  usually 
in  a  secular  manner j .  Sea  Reference  20. 

2.  Fluctuations  in  the  aolar  constant  produce  stoehaatlc  varia¬ 
tions  about  the  average  solar  pressure  force  described  just 
above.  These  fluctuations  are  not  very  predictable  and 
anount  to  perhaps  .1  to  .2 i.  On  a  Jisriner  class  spacecraft 
(01.5  lbs.)  this  would  produce  acceleration  variations  on  the 
order  of  10*1’  kre/se c2  at  Hers  dlster.ee.  See  Reference  20. 


3.  Oas  leakage  from  the  *eta  on  the  attitude  control  system  have 
proved  the  most  troublesome  In  the  post.  The  leakage  is  due 
to  such  phenomena  aa  dust  particles  trapped  in  the  valves 
themselves  (this  changes  each  tine  the  valves  fire  which  la 
perhaps  once/hr  and  in  general  tends  to  be  erratic)  and  im¬ 
balance  in  the  couplinp.  Each  spacecraft  flown  shows  differ¬ 
ent  characteristics  in  this  area.  To  modal  it  requires  care¬ 
ful  study  of  the  attitude  control  sub-system  daU  and  is  not 
really  tracteble,  particularly  during  the  real-tine  mission. 
See  Reference  21. 


*».  Leakage  fror  the  propellant  systems  and  scientific  instru¬ 
ments.  This  item  is  highly  epacecrrft  dependent  and  will  be¬ 
came  more  important  aa  spacecraft  carry  more  propellant  under 
high  pressure  (c.g.,  for  planetary  orbitera)  oni  sore  caa- 
plex  scientific  packages  (e.g.,  the  infrared  spectrometer 
gas  blowdown  during  Mariner  V2  end  VIZ  encounter  blurred  the 
information  inherent  in  the  tracking  dots  enough  to  cake  it 
nearly  unusable  for  celestial  aechanlcs  purposes).  8ee 
Reference  22. 


iok 


3.  For  solar  electric  or  nuclear  propula  ion  ijrtUae,  the  powered 
flight  phase  extends  during  the  entire  (or  nearly)  cmiee 
phaae.  Fluctuations  in  the  thrust  Magnitude  and  direction 
produce  stochastic  spacecraft  accelerations  perhaps  as  high 
as  10*°  ta/acr  .  These  are  so  high  as  to  asks  the  analysis 
of  the  situations  qualitatively  (as  well  as  quantitatively) 
different  fron  the  so-called  ballistic  cases  we  are  studying 
and  are  beyond  the  scope  of  this  effort.  See  Jordan 
(he  Terence  23) 

The  preceding  list  was  for  background  and  to  convince  that  al¬ 
though  each  specific  problaa  had  different  isplicatioos  and  should  be 
studied  with  a  specific  a  lesion  and  spacecraft  in  Bind,  cannon  to  all 
was  that,  whatever  the  situation,  our  current  knowledge  does  not  pro¬ 
vide  an  accurate  stochastic  nodal  for  the  spacecraft  accelerations 
with  which  to  derive  o  aUlasni  variance  filter,  tone  structure  la 
needed  and  n  nodal  thoegfct  lihsly  for  n  Viking  '73  spacecraft  will 


he  aseunsdi 


1.  There  will  he  e  fairly  larva  unknown  bwt  constant  force 
arising  fron  the  leakage  in  the  kigh  preeeure  propellant  storage  areas. 
This  will  be  tabes  a-  a  naan  aero,  spherically  distributed  accelera¬ 
tion  cf  standard  deviation  10*U  lessee2  in  oooh  of  the  tkrue  ortho 


f .  In  addition,  there  will  be  the  rnadea  forces  fron  ths  neurone 
Juet  described,  for  the  bnselln  nodal  need  bare,  a  standard  devia¬ 
tion  of  10’u  ha/eoc2  will  bo  aaeunsd  and,  again,  a  correlation  tine 
of  29  dope,  jsrsnss  the  randan  colons  ate  are  not  eetiaated,  the  nag- 
ml  tads  af  ths  acosleratloa  sasid  is  not  lap  onset  -  i.s.,  an  aceel- 
eration  af  trios  tbs  standard  dsviatisn,  for  esaaple,  will  produce 
twice  the  effect.  The  correlation  tint  is  Ufortaat  and  supporting 
studios  art  psrfosasd  whore  this  is  varied. 
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»>mu  -  C»»»  1  -  and  long  Arc  Powlf  <*\1 


The  a  mm  Viking  trajectory  u  mi  m«4  for  the  a  tat  l  on  location 
•tody  mm  selected.  figures  6*6  and  6.7  show  bene vl or  ortr  a  50-day 

a 

period.  Displayed  la  tlM  uncertainty  added  la  t be  r  direction  by  tba 
presence  of  tha  atocbaatlc  forcta.  An  independent  process  it  aasuned 
for  each  of  tha  three  orthogonal  spacecraft  axes.  Of  course,  fron 
tba  previous  analysis,  only  tba  eonponant  of  acca  la  ration  in  tba 
larth  probe  direction  la  expected  to  bare  nach  affect  and.  In  the  ab- 
aanea  of  ranging  data,  only  tha  radial  ecnponant  of  tba  poaltloo 
oatlnnta  would  be  largely  affected. 

Tba  batten  curve  In  figure  6.6  la  for  a  correlation  tine  tf  ID* 
daya,  or  oaaantlally  a  oonatant  proeeaa.  Baaed  on  tba  meal  to  of 
OMftor  X¥,  a  oontrlbntlon  to  tba  aneartalaty  In  r  of  abont  50  kn  la 
anpaetod  (we  are  aalng  10*12  In^aoc2  lnataod  of  10*U  bn/aoc2  even 
fbr  tba  oonatant  forao  to  faalUtato  caagerlaon).  But  Jvnt  aa  In  tba 
atatlon  loaatlon  eaao,  tba  analytic  conpntatloa  la  not  accurate  for 
Zt  rlaoo  gal  JOy  Iran  tba  50  bn  value,  peaking  at  200  ha  after 
•bant  10  daya  of  tracking.  Tha  analytic  foraala  la  a  factor  of  h  la 

Bean  wore  disturbing  la  that  tba  errora  Wild  rapidly  aa  tha 
correlation  tint  la  docreaaad.  Tba  blghact  contribution  observed  la 
for  a  1-day  correlation  tint  after  >5  daya  of  tracking  •  1500  kn. 

Tbla  la  a  factor  of  50  higher  than  war  calculated  analytically  for  tba 
conataat  acceleration  caae.  fortunately,  aa  sore  data  la  taken,  the 

affect  dlnlnlahee  until  after  20  daya  there  la  little  difference  be¬ 
tween  ▼  •  1  end  t  ■  25*  In  c*neral,  t  »  5  appears  to  be  the  worst 
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HACKING  TIME,  doyi 


Figure  6.6.  Effect  of  Random  Acceleration*  on  8Ute  Only  Filter 

venue  Tracking  Arc,  T,  For  Several  Different  Correlations 
Doppler  Only;  Lx  is  Earth-Probe  Position  Component. 


HACKING  TIME,  do y* 

Fifure  6.7.  Effect  of  Rendon  Accelerations  on  State  Plus  Three 
Constant  Accelerations  Filter 
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correlation  time  possible. 

Figure  6.7  plots  ths  saae  information  for  s  filter  that  estimates 
three  constant,  orthogonal  accelerations  assuaing  an  aprlorl  of  10* 11 
km/sec.  The  perforaance  is  not  auch  better;  in  fact,  after  about  20 
days  or  so,  all  processes  plotted  have  greater  effect  on  the  extended 
filter  than  they  do  on  the  restricted,  state-only  filter.  This  is  not 
to  say,  however,  that  the  overall  perforaance  is  not  laproved  since 
the  extended  filter  handles  any  conatant  forces  that  are  present  quite 
wall.  Where  there  actually  are  constant  accelerations  present  at  the 
10  ka/sec  uncertainty  level  and  10  ka/sec  stochastic  accelera¬ 
tions,  the  overall  perforaance  (not  shown)  was  better  in  every  case 
studied. 

L 

Figure  6.8  shows  the  effect  of  the  t  ■  10  and  5- day  accelera¬ 
tions  on  a  state  only  filter  over  a  auch  longer  six-aonth  tracking 
Interval.  The  trends  are  soswidiat  erratic  but  the  aain  aessage  is 
clear-process  noise  over  the  data  arcs  typical  for  inner  planet 
■lesions  does  not  produce  filter  divergence. 


Results  -  Case  2  -  Medium  Arc  with  and  without  Range  Data 

The  preceding  showed  a  remarkable  sensitivity  to  stochastic 

forces.  Although  we  consider  it  quite  illuminating,  the  cases  studied 

1 

contain  two  difficulties: 

1.  No  range  data  were  simulated  -  it  was  found  in  Chapter  IV 
that  range  data  were  quite  useful  in  reducing  the  unmodeled 
force  effects,  particularly  when  the  estimate  did  not  have 
to  be  propogated  very  far. 

2.  To  keep  the  amount  of  data  presented  manageable,  only  the 
rango  component  of  position  has  been  discussed.  With  rang¬ 
ing  data,  velocity  as  well  as  position  errors  are  induced, 

but  presenting  results  in  all  fl  'c  coordinates  involved  is 
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T,  doy* 

Figure  6.8.  long  Term  Effect  of  Random  Accelerations  on  State 
Only  Filter 
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•  very  awkward* 

To  circumvent  both  difficulties,  we  shall  now  consider  a  case 
including  ranging  data  and  display  the  results  in  tens  of  the  norm 
of  position  of  the  propogated  state  error.  Specifically,  a  new 
Viking  '7 5  trajectory  will  be  considered  with  epoch  and  data  starting 
at  Kars  encounter  minus  30  days  (E-30d) .  Tracking  data  is  simulated 
continuously  from  three  stations  until  E-10  days  at  which  time  the 
best  estimate  possible  is  required  in  order  to  calculate  a  last  tra¬ 
jectory  correction.  Since  the  mission  criteria  will  require  delivery 
of  the  spacecraft  to  a  certain  Mars  relative  position,  the  uncertain¬ 
ties  contributed  by  the  random  accelerations  are  prorogated  to  encoun¬ 
ter  and  the  square-root  of  the  trace  of  the  position  covariance  is 
calculated  at  that  time.  In  this  manner  a  realistic  scalar  measure  to 
the  6  dimensional  filter  performance  can  conveniently  be  assigned. 

This  performance  measure  is  shown  plotted  in  Figure  6.9  for  four 
separate  cases  which  Include  the  combinations  of  ranging  data  and  not, 
and  extended  three  constant  force  solution  and  not.  Each  case  assumes 
a  correlation  time  on  the  stochastic  processes  of  10  days.  Note  that 

the  performance  measure  is  not  the  overall  filter  performance,  it  is 

,  / 

the  effect  of  the  stochastic  accelerations  only.  The  following  obser¬ 
vations  may  be  made: 

1.  Using  doppler  data  only,  the  contribution  from  the  forces 
assumed  is  roughly  1000  km  whether  an  extended  solution  is 
made  or  not. 

2.  The  addition  of  ranging  data  to  the  state  only  case  improves 
the  performance  for  very  short  arcs,  but  by  the  time  20  days 
of  data  hove  been  accumulated,  the  error  again  has  risen  to 
approximately  1000  km. 


Figure  6.9.  Norn  of  RMS  Position  Error  at  Encounter  Resulting 
from  Stochastic  Accelerations. 
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3.  Using  both  ranging  data  and  extending  the  solution  produces 
a  markedly  superior  result.  After  20  days  of  tracking,  the 
error  induced  is  only  170  km. 

It  is  veil  to  check  the  performance  of  this  later  case  to  make 
sure  that  the  results*  are  not  highly  sensitive  to  the  assuuwd  corre¬ 
lation  time  of  10  days.  Figure  6.10  shows  results  for 

1 

5 

t  ■  10  days 

25 

using  ranging  data  and  the  extended  solution.  As  Bight  be  expected, 
performance  improves  for  the  longer  correlation  times  since  the  pro¬ 
cess  is  more  nearly  a  constant  and  hence  can  be  removed  efficiently  by 
the  extended  pnru:.ieters  in  the  solution.  The  exception  to  this  is  the 
T  ■  1  day  case  which  after  20  days,  produces  less  effect  than  the 
▼  ■  5  day  case.  This  is  to  be  expected  because  in  the  Halt  of  zero 
correlation  time,  the  total  power  o**  the  stochastic  process  will  be 
zero.  Thus  the  situation  arises  where  there  should  be  a  "wo’-st" 
correlation  time  and  the  results  Indicate  we  have  found  it  -  at 
t  M  5  days. 

Finally,  Figure  6.11  presents  data  yielding  the  total  filter  per¬ 
formance  including  the  effects  of  data  noise  and  the  constant  forces. 
There  arc  three  curves  that  can  meaning!'  lly  be  colculatcd: 

©  Effects  of  data  noise  alone,  i.c.,  the  filter  solves  for 
three  constant  forces  but  there  are  actually  none  present. 

(§)  Data  noise  plus  constant  accelerations,  i.e.,  the  filter 
is  the  srnc  but  this  time  the  constant  accelerations 
actually,  exist  and  at  the  level  assumed  by  the  filter 
(o  -  10"^  km/sec^) 
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Figure  6.10.  Worm  of  RHS  Position  Error  at  Final  Tirre 
Resulting  from  Stochastic  Accelerations. 


©  Data  noise,  constant  accelerations,  and  stochastic 
accelerations. 

It  is  scon  that  the  three  curves  are  nearly  indistinguishable 
for  tracking  arcs  of  less  than  eight  days.  After  that  point,  they 
separate  until,  at  the  end  of  the  tracking  arc,  curve @  is  more  than 
40 %  higher  than  curve Q.  The  A  is  the  final  performance  of  an  idea¬ 
lized  case  where  no  accelerations  of  any  kind  are  present,  permitting 
a  state  only  solution. 

Some  General  Considerations  of  Worst  Cose  Analysis 


Random,  unmodeled  process  noise  has  been  shown  to  have  large 
effects  on  state  estimation  accuracy  even  when  present  at  the  extrem¬ 
ely  small  levels  assumed.  In  the  course  of  this  Investigation,  we 
have  resisted  believing  the  results  (because  the  effects  are  so  pro¬ 
nounced)  and  have  developed,  we  believe,  illuminating  points  of  view 
in  attempts  to  understand  better  the  origin  of  the  results  obtained. 

In  Chapter  III  it  was  derived  that,  for  a  few  days  at  least,  the 
partial  of  the  doppler  with  respect  to  geocentric  range  could  be 


written 


iL 

a* 


kt 


Since  a  constant  acceleration  would  affect  the  data  as 


AP  -  k't 

it  would  be  indistinguishable  from  r.  This  is  not  to  say,  that  a  con¬ 
stant  acceleration  would  produce  the  max  inrum  error  in  r,  however.  In 
order  to  see  this,  suppose  that  a  variable,  x^,  is  being  estimated 
whose  data  partial  over  the  available  data  set  is  orthogonal  to  every 
other  data  partial  of  the  remaining  variables  being  estimated.  This 


. 


■H 
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partial  Is  shown  in  a  two-dimensional  sub-space  of  the  m  dimensional 
space  spanned  by  the  m  data  points  being  taken  (Figure  6.12).  In  the 
Householder  notation  ured,  R  will  have  the  form 


Figure  6.12.  Maximum  Error  Occurs  when  n  lies  along  r^. 

In  this  circumstance,  the  error  in  v*U  be  simply  the  dot 
product  of  n  on  the  x1  partial.  If,  for  example,  the  worst  possible 
noise  of  a  given  norm  was  sought,  it  would  dbviously  be  that  n  that 
was  co-llnear  with  the  original  date  partial. 

Suppose,  however,  that  an  additional  parameter  wtue  added  to  the 
estimate  list,  vhose  partial  was  contained  within  the  e^,  sub* 
space  just  described.  Coll  that  partial  rR4^  -  it  is  shown  in 
Figure  6.13. 
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Figure  6.13.  When  the  data  partial*  ere  not  orthogonal,  the 
Maximum  error  does  not  occur  when  n  lies  along 

V 

How  the  portion  of  the  noise  that  lies  in  this  sub-space  must  be 
allocated  to  the  two  estimates  6^  and  $n+1;  i.e.,  none  of  the  noise 
will  be  placed  in  the  residuals.  For  exaaple,  if  n  lies  along  r  .  , 
no  error  will  be  induced  in  x^.  The  Maximus  error  in  is  induced 
when  n  is  orthogonal  to  This  circumstance  is  shown  in 

Figure  6.lU. 


ore  highly  correlated.  Going  bach  to  the  analysis  of  Hunter  TIT, 


•Mfut 

t 

with  tin  t 
and  t  co*  t 

These  functions,  and  hence  the  data  partial*  art  vary  siailar  - 
h*nce  the  aaxlaua  error  will  not  be  caused  by  a  noise  vector  in  the 
direction  of  the  original  partial. 

Consider  a  second  point  of  viev.  Suppose  a  standard  least 
squares  estimator  is  employed: 

*  .(*VlA 


and  suppose  that  the  A  matrix  ia  given  by 


now  if 

a  ■  Ax  4  n 
and  definite 

r  ft  (ATA)‘lAT 

then 

ft  -  x  •  Fn 

The  F  Matrix  is  readily  seen  to  be 


Suppose  the  error  in  ft^,  only  is  of  direct  interest.  It  is  desired 
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to  calculate  the  largest  error  In  this  variable  if,  say,  each  element 
of  n  were  bounded  between  +1,  but  was  arbitrary  otherwise.  If 


were  selected  (i.e.,  a  noise  aignature  that  "looked"  like  the  original 
data  partial),  tho  error  would  be  unity.  If  however 


were  selected,  en  error  of  1.9  would  result,  the  largest  possible 
for  the  problee  posed. 

To  return  to  the  process  noise  situation  under  study,  suppose 
the  vorst  possible  uncodcled  acceleration  signature  is  desired  with 


the  only  constraint  being  that  the  acceleration  be  a  bounded  variable. 
Consider  the  row  of  the  K  matrix  of  the  coordinate  being  studied  (r 
would  be  ideal)  as  a  continuous  function.  The  error'  Induced  would  be 


where  is  the  perturbation  in  the  data, 
(assuming  doppler  data  only)  as  (t) 


If  the  approximation  that 


a  (u)  du 


is  awde  where  a  (u)  is  the  bounded  acceleration  (this  would  be  en 
excellent  spproxiastlon  on  a  deep  space  trajectory  for  a  month  or  so). 


The  problem  becor.es: 
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Extreaalize  subject  to 

|  •  (t)  |  <  K 

This  is  a  straightforward  Calculus  of  Variation*  problaa.  Define  tin 
Haailtonian 


H  (**,  a,  l»  t)  -  f  (t)  ftf.  (t)  ♦  x(t)a(t) 
Then  Mlniairlnc  with  respect  to  a  (t),  we  find 
a  (t)  -  -K  sen  \  (t) 

8incc  on  an  cxlrcnal  arc 


Than 


-fe- 


»(t) 


H 


8* 


-  -  i(t) 


■•r 


f(r)  dt  ♦  C 


Bacauaa  there  are  no  end-point  constraint* 
X  (t)  -  0 


and 


thus 


X  (t) 


/t 

J  f  (r)  dr 
/ 


•  (t)  •  -K  *gn  (  f  (▼)  dt  *  (i) 

« 

This  is  the  answer  in  the  continuous  fomuiation  -  the  discrete 
version  yields  highly  siailar  answers.  Figure  6.15  plots  the  discrete 


«tfo  are  following  n  straightforward  application  of  Fontryagin's 
Kcxiaua  principle  as  outlined  in  D.  H.  Uiberg's  class  notes  for 
Engr.  2??-C,  Winter  Quarter,  1969.  The  development  is  part  of  a 
separate  project  being  carried  out  by  Jataeu  Mcbonell  and  this  author. 
It  is  beyond  the  scope  of  this  thesis,  but  we  wish  to  apply  some  of 
the  results  here  and  will  deacrib*  the  procedures  being  used  briefly. 


* 


*>1 
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Figure  6.13.  Discrete  F  Function  with  respect  to  r  coordinate  for  eight  days  Doggier  Trechl 


f  functions  for  the  r  estlMtc  for  tht  first  eight  days  of  dnppler 
data  on  the  Viking  trajectory  of  Csss  Z  Just  abort.  Kota  how  dia- 
•inllar  it  is  fr on  the  single  k*t  structure  of  the  original  partial. 
If  equation  (I)  is  executed  by  sunning  the  points  of  this  plot  start* 
Ing  at  the  rlfM-htrH  side,  there  art  two  reversals  of  sign  as  the 
manlag  is  accusal e ted.  The,  the  worst  com  acceleration  Is  one 
thst  starts  at  ♦*,  reverses  once  to  -K,  and  finally  reverts  to  «K 
scale  es  T  interval  la  t reversed.  This  is  hardly  a  constant  accelera¬ 
tion.  Thin  example  lends  credence  to  the  covariance  results  tdtich 
determined  that  e  5-day  correlation  tine  is  the  worst  case  end  given 
•one  insight  as  to  the  origin  of  thie  result. 
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CHAITER  VII 


FURTHER  CONSIDERATIONS  IN  DETECTING  UNMODELED  PARAMETERS 

Much  of  the  preceding  material  has  dealt  with  filter  response  to 
unaodeled  or  non-precisely  modeled  parameters  and  processes.  The 
presence  of  un-modeled  parameters  can  seriously  degrade  state  estimates 
as  has  been  shown.  Solving  for  them,  even  non-optimally,  can  often 
be  very  helpful  as  was  detailed  in  Chapter  VI.  There  is  a  price  for 
this  solution,  though,  and  that  is  that  the  state  estimates  are  more 
sensitive  to  data  noise.  That  is,  as  more  parameters  are  included  in 
the  solution,  the  filter  becomes  more  broadly  tuned,  and  its  discrim¬ 
ination  powers  between  signal  and  noise  are  correspondingly  less. 

Thus,  if  a  parameter  is  thought  not  to  be  present  (at  a  significant 
level)  but  it  is  feared  that  it  night  be,  the  "fear"  can  be  covered 
(often)  with  a  solution  but  given  that  the  "thought" was  correct,  the 
covariance  of  state  errors  will  be  higher  if  a  solution  is  made  than 
if  it  is  not. 

These  remarks  can  be  put  on  a  quantitative  basis  by  employing 
the  generalized  consider  option  of  Chapter  V.  For  example,  filter 
performance  can  be  calculated  for  the  situations  where  constant  non- 
gravltatlonal  forces  are  either  solved-for  or  not  and  are  present  or 
not.  Table  7.1  displays  the  results  for  4  separate  cases  for  two 
different  epans  of  doppler  tracking  data: 

1.  Filter  estimates  state  only;  forces  not  present. 

2.  Filter  estimates  stete  and  three  forces;  but  forces  are  not 

Preceding  page  link 
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Physical  Model  and  Filter  Design. 


present. 

3.  Filter  estimates  state  and  three  forces;  forces  do  exist. 

4.  Filter  estimates  state  only;  forces  exist. 

The  variable  recorded  is  the  r  direction,  the  main  state  parameter 
affected  by  unmodeled  forces.  In  general,  while  the  true  filter  per¬ 
formance  follows  the  numbering  system  Just  adopted,  the  specific 
point  being  made  is  that  Case  2  shows  degraded  performance  from  Case  1. 

The  two  parameters  classes  dealt  with  in  previous  chapters 
(station  locations  and  non-gravitational  forcer;)  present  the  very  dif¬ 
ficulties  for  real  time  orbit  determination  of  which  we  are  speaking. 
Spacecraft  are  designed  carefully  to  minimize  translational  forces 
being  spuriously  generated  by  the  attitude  control,  propulsion,  and 
other  systems  for  this  very  reason.  They  can  develop  (and  have,  see 

O,*  lii1  *.  ,  x  »  1  •*  •*  ■*»  .  j  '  '  i  » * ,  i  .  ( ■*  *  *  >.  j’ 

Ref. 26)  leaks  at  any  time,  however.  As  for  station  locations,  these 

**  J  f-  i  * 

parameters  are  carefully  determined  using  data  from  previous  missions- 
there  is,  however,  the  not  very  remote  possibility  that  the  solutions 
are  not  as  accurate  as  is  thought,  the  propogation  model  is  inaccurate, 
or  that  the  universal  time  calibrations  are  inaccurate  for  the  new 
time  spans  being  used. 

These  factors  combine  to  motivate  the  search  for  computational 

machinery  to  augment  existing  orbit  determination  methods  that  aid  the 

-  .  j  .  >■ 

analyst  in  performing  formal  or  informal  hypothesis  testing  with  re¬ 
gard  to  his  solutions.  The  data  which  he  has  at  his  disposal  are 
l)  the  temporal  or  data  set  dependent  behavior  of  the  solutions  them¬ 
selves  and  2)  the  data  residuals  after  the  solution(s)  has  been  made. 
In  this  chapter  we  wlnh  to  surest  a  simple  algorithm  to  serve  as  a 


"detectability  index",  a  means  of  deciding  to  what  extent  an  unmodeled 
parameter,  if  present,  would  be  noticed  by  Inspecting  the  data  resid¬ 
uals,  or  more  simply,  the  weighted  sum-of-squares  of  the  residuals. 

A  Considered  Residuals  Option 

Calculating  the  signature  impressed  on  the  data  by  a  hypothetical 
parameter  is  not  directly  useful  in  pursuing  the  goal  of  detecting 
unmodeled  parameters.  To  look  for  this  signature,  one  would  have  to 
inspect  the  data  before  the  fit,  or  solution,  were  made  and  this  data 
would  contain  the  signatures  from  all  other  state  parameters  as  well. 
We  must  be  content  looking  for  this  desired  signature  after  it  has 
passed  through  the  filter  to  be  used.  Of  course  its  character  will 
be  changed  by  this  operation  -  formally,  what  will  remain  is  the  por¬ 
tion  orthogonal  to  the  column  space  of  the  partial  matrix  used  to  de¬ 
sign  the  filter. 

For  example,  consider  a  data  set  linearly  related  to  the  estimated 
state,  x,  and  unmodeled  parameters,  y,  as  follows 

z  ■  Hx  +  By  +  c  (7.1) 

Clearly,  y  adds  By  to  z.  If  the  least  squares  filter, 
it  -  (H1  A  “by1  H1  a/1  z  -  F  z 

is  employed  and  !&  is  subtracted  from  z  to  form  the  residual  vector, 
then 

«  :  '  '  0  '  1  •  '  »  *  /  '  * 

*z  -  z  -  Hk  -  (I  -  H  F)  By  +  (I  -  H  F)  a  (7.2) 

The  first  half  of  the  right  hand  portion  of  7.2  are  the  residuals  due 
to  the  presence  of  y. 

In  any  given  problem,  having  on  hand  this  filtered  signature  to 
compare  with  the  residual  patterns  actually  obtained  might  be  quite 
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useful.  It  has  the  disadvantage  of  being  unwieldy,  unfortunately, 

. 

since  its  application  would  require  plotting  of  these  signatures  for 
each  parameter  being  "feared"  and  for  each  data  set  (or  span)  used  for 
the  various  solutions  made.  The  curse  of  too  much  information  lurks 
nearby  this  suggestion.  Something  of  more  facile  application  would  be 
to  calculate  the  increase  in  the  weighted  sum-of-squeres  due  to  the 
existence  of  an  unmodelcd  parameter.  Using  7.2  this  is  easily  seen  to 
be 

8zTgz1  due  toy  -  yTBT  (I-H  F)T  (I-H  F)  By 

or  (7.3) 

«  trace  (I-H  F)  By  yTBT  (I-H  F)T 

Unfortunately,  either  form  o**  (7.3)  involves  literal  manipulations 
with  matrices  of  very  inconvenient  dimensions  (for  example  the  number 
of  rows  of  B  is  equal  to  the  dimension  of  the  data  vector).  The  House¬ 
holder  transformation  technique,  which  automatically  selects  a  new 
basis  for  the  partial  derivative  matrix  can  be  conveniently  employed 
to  solve  the  dimensionality  problem  imbedded  in  7. 3* 

We  begin  the  formal  development  by  once  again  employing  7.1  re¬ 
cast  as  a  square-root  filtering  problem  as  described  in  the  Appendix 
and  extended  in  Chapter  V. 


As  before,  a  solution  for  7.1*  is  desired  that  minimizes  the  norm  of 
the  weighted  residual  vector.  The  y  variables  are  again  the  unmodelcd 
parameters  which  may  or  may  not  be  present  -  their  estimates  will  be 
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left  at  zero  in  either  case.  Pre-multiplying  7.^  by  the  square-root, 
weighting  matrix  and  then  performing  the  Householder  transformation, 
recall  the  results  of  Chapter  V  and  re-write  as: 


(7.5) 


Note  that  the  presence  of  y  results  in  a  contribution  to  the  residuals 
only  in  p  of  the  elements  in  terms  of  the  new  basis  which  has  been 
created  (where  p  is  the  dimension  of  y).  It  is  convenient  to  think 
of  y  os  a  zero  mean  random  variable  uncorrelated  with  x  (as  was  done 
in  the  development  of  the  consider  option)  and  compute  the  covariance 
of  £-x,  fty,  Post-multiplying  7.5  by  its  transpose  and  taking 
ensemble  averages  . 
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(7.6) 


The  upper  left  partition  of  7.6  is  the  covariance  of  state  es¬ 
timate  errors  and  is  identical  to  5.3.  Writing  out  the  lower  right 
2x2  partition,  the  covariance  of  residuals: 


VA*  + 1  0  }  p  (7.7) 

0  I  }  m-(p+n) 

where  p  it  the  dimension  of  the  consider  parameter  vector 
m  is  the  dimension  of  the  data  vector 
n  is  the  dimension  of  the  consider  solution  vector 
and  since 

T  T 

xx"  trace  x  x 

it  is  concluded  that  the  trace  of  (7.7)  is  the  expected  value  of  the 
weighted  sum-of-squares  of  the  residual  vector.  Note  that  if  A  is 

y 

zero,  this  is  simply  m-n,  the  classical  least-squares  result  for  uni¬ 
formly  weighted  data.  It  is  interesting  that  this  is  also  the  result 
for  an  arbitrary  weighting  of  the  data,  W,  as  was  employed  here.  More 
compactly: 

E^ weighted  sum-of-8 qua real  A  sos  -  (m-n)  +  tr  RyAyRyT  (7.8) 

Equation  (7.8)  can  be  computed  efficiently  and  is  a  useful  result  in 
Itself  but  still  does  not  get  at  the  notion  of  detectability.  It 
would  be  more  desire  able  to  compute,  say,  the  per-cent  increase  in  sos, 
per  unit  change  in  the  estimate  due  to  an  unmodeled  parameter.  Equa¬ 
tion  (7*8)  is  cast  statistically  but  it  is  obviously  possible  to  write 
the  Increase  in  the  sos  in  terms  of  a  particular  y,  i.e., 

AS  os  -  tr  Ry  y  yT  RyT 

For  example,  if  p  were  2  then 

Asos  -  ru2  yx2  ♦  (r^2  +  r^2)  yg2  +  2^^  y^  (7.9) 
where  r^  is  the  i,  J  element  of  Ry,  an  upper  triangular  matrix. 
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Generalizing  (7.9)  it  is  easily  shown  that  the  Increase  in  the  sos  due 
to  the  presence  of  a  single  parameter  y^  is 

"°*l  *  }i irJi2  (7>l0) 

The  cross-product  terms  do  not  appear  in  (7*10)  since  each  parameter 
is  treated  separately  -  that  is,  in  postulating  that  y^  has  a  non- zero 
magnitude,  it  is  implicitly  assumed  that  the  rest  of  the  parameters 
are  zero. 

As  for  the  sensitivity  of  the  estimate  to  the  unmodeled  parameter, 

«  • 

this  ia  readily  seen  to  be  (from  7.5) 

*  ’  V1  V 

Again  considering  only  one  element  of  y,  y^, 

A*  -  Rx’1  Txy  (i)  yA  (7.11) 

where  (i)  ia  the  ith  column  of  R^. 

Finally,  by  combining  (7.10)  and  (7.11),  the  following  suggested 
inverse  detectability  index  is  established. 

di  -  V1  rxy  “>/  VA 

In  words,  (7.12)  gives  the  vector  change  in  the  estimated  para¬ 
meters  divided  by  the  per  cent  increase  in  the  rss  weighted  residuals 
(over  that  expected  statistically  from  noise  alone)  due  to  the  ith 
unmodeled  parameter.  Thus,  speaking  qualitatively,  in  a  real  fitting 
problem,  even  though  it  can  be  demonstrated  that  the  estimate  is  sen¬ 
sitive  to  y^,  if  d^  is  small,  and  the  residuals  conform  well  to  that 
expected  from  noise  alone,  it  can  be  concluded  at  c  glance  that 
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yi  is  not  present  at  a  damaging  level. 

This  index  should  work  veil  in  validating  apparently  good  solu¬ 
tions  in  a  real-time  situation.  If  the  residuals  do  not  show  struc¬ 
ture,  the  sos  is  nearly  what  is  expected,  and  all  the  d^'s  are  adequat¬ 
ely  low,  then  it  could  be  concluded  with  reasonable  assurance  that 
nothing  untoward  has  occurred  and  the  modeling  has  been  faithfully 
executed.  In  pre-cnalysis,  the  index  might  be  useful  in  establishing 
what  data  sets  are  appropriate  for  the  execution  of  the  validation 
process.  In  general,  the  d^fs  should  become  lower  with  increased 
amounts  of  data. 

If,  on  the  other  hand,  the  date  does  not  fit  well  and  there  are 
several  parameters  that  could  be  causing  the  difficulty,  it  would  be 
useful  to  have  a  convenient  means  of  investigating  which  of  these 
parameters  (if  any)  when  included  in  the  solution  would  decrease  the 
sos  by  significant  amounts.  It  would  be  desireable  to  have  an  algor¬ 
ithm  to  do  this  without  explicitly  performing  the  various  solutions*. 
To  find  the  desired  algorithm,  first  break  the  consider  parameter 
vector,  y,  into  three  sub-classes: 


In  terms  of  these  newly  defined  y^'s,  we  wish  to  compute  the  reduction 

on  the  scs  that  would  be  obtained  if  y?  were  added  to  the  estimation 

•That  this  could  be  done  was  first  pointed  out  to  the  author  by  S.  R. 
McReynolds  in  a  private  communication. 
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list,  but  the  estimate*  of  and  were  still  left  at  zero.  Again, 
yg  is  restricted  to  having  unity  dimension  and  only  the  possibility  of 
adding  a  a ingle  new  variable  to  the  estimation  list  is  investigated. 
The  following  equation  must  be  solved: 
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(7.13) 


for  x,  ^  *uch  that  the  norm  of  A  is  minimized. 
Writing  (7.13)  as  a  set  of  equations: 


R  x  +  r 


xy. 


\y2h  *  \ 


(7.1U) 


f2  *  y2 


From  the  structure  of  7.1k: 

1.  There  is  no  freedom  to  alter  n^  fro n  that  specified. 

a  a  3 

2.  Whatever  the  choice  of  /g,  *  car.  be  adjusted  so  that  nx 
may  be  zero. 

3.  Therefore  the  problem  resolves  to  one  of  choosing  y?  so 


13<l 


that  th« 


la  a  minimum 
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From  (7.ll4),  it  la  evident  that 


ny^T  ny^  ♦  ny^2  -  (zy^T  zy^  ♦  *y^2)-2(zy^T  r  +  ^  ry^)yg 
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Minimizing  (7*1?)  with  reapect  to  y_  yielda 


(7.15) 


(*y  T  ry  y 

yl  yly2 


♦  z  r  ) 
yl  y2 


(r  T  r 
yly2  yly2 


*r  2) 


and  the  decrease  in  the  a  os  is  simply 


AS  os 


yl  yly2  yl  y2 

«  p 

(p  r  ♦  p 

yly2  yly2  y2 


(7.16) 


As  an  application  example,  suppose  that  several  stations  are 
taking  data  and  through  a  blunder  a  single  component  of  one  of  the 
station  locations  is  seriously  in  error  and  none  of  these  locations 
are  Included  in  the  solution.  All  of  the  data  residuals  will  be 
affected  by  the  error  as  a  best  compromise  is  found  by  the  lease* 
square  algorithm.  The  actual  error  component  could  be  Isolated  very 
quickly,  however,  by  executing  (7.1iS)  for  each  consider  peraueter, 
since  only  that  paraaeter  would  be  completely  successful  in  removing 
the  structured  residual  patterns  that  have  been  produced. 


A  Numerical  Example 


An  1106  Computer  program  vai  mechanized  to  calculate  the  a  os 
changes  and  the  detectability  index  for  several  possible  modeling 
errors.  Again  as  in  the  preceding  chapter,  the  program  was  designed 
to  use  the  partials  generated  by  the  DPODP.  To  continue  the  theme  of 
the  earlier  chapters,  non-gravitational  forces  and  the  station  loca¬ 
tion  errors  are  selected  for  the  numerical  examples. 

Figure  7*1  plots  the  increase  in  the  sos  over  that  expected,  from 


data  noise  alone,  l.e., 


EfsosJ 


E  r2/(m-n) 


assuming  a  solution  for  the  six  components  of  the  trajectory  was  made. 

The  DFODP  models  the  non-gravitational  forces  as  components,  not 
in  the  Earth- relative  plane-of-the-sky,  but  along  the  spacecraft  axes. 
Since  the  roll  axis  is  oriented  along  the  sun-probe  line,  this  be¬ 
comes  one  component,  ar.  The  other  two  componets  are  ay  which  is  in 
the  ecliptic  but  perpendicular  to  the  direction  of  ar  and  a  third 
direction  4hlch  is  perpendicular  to  the  ecliptic,  ax.  The  results 
could  be  related  better  to  the  preceding  chapters  if  the  plane-of-the- 
sky  coordinate  system  were  used,  but  we  are  interested  in  detecting 
general  classes  of  parameters  hcre-not  Just  those  which  are  well 
understood.  This  rather  arbitrary  coordinate  system  will  suffice. 

The  station  location  errors  are  denoted  in  the  coordinate  system 
being  used,  rg  and  X. 

From  Figure  7.1,  as  expected  the  per  cent  change  in  the  sos  rises 
as  more  trsckinr  datr  are  included.  The  .forral  principle  at  work  here 


136 


ASOS/E(SOS) 
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Figure  7.1.  Change  in  sos  divided  by  that  expected 
versus  Tracking  Arc. 
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is  that  with  additional  data,  the  data  set  spans  more  completely  the 
consider  parameter  space  as  well  os  the  solution  space  and  data  caused 
by  the  unmodeled  parameters  would  be  more  efficiently  forced  out  of 
the  solution  and  into  the  residuals.  Of  the  three  accelerations,  a 

w 

is  by  far  the  most  visible  in  the  data  (the  y  direction  most  nearly 

r 

coincides  with  the  Earth-line-of-sight) .  For  very  short  spans,  an 
rg  error  is  more  visible  than  a  \  error.  This  relationship  reverses 
after  a  few  days  until,  after  30  days  of  tracking,  a  10m  longitude 

fit;  *r  , 

error  will  increase  the  residuals  by  Uo£,  whereas  the  same  magnitude 

*  *  i  •  . 

of  r  displacement  yields  but  a  7 %  increase.  The  speculation  is  made 
s 

that  this  is  because  the  right  ascension  of  the  spacecraft  (which 
appears  in  the  data  like  a  station  longitude)  is  influenced  more  by 
the  solar  gravitation  than  is  declination,  creating  a  greater  separa¬ 
tion  between  >.  and  u  than  between  rg  and  5. 

Equation  7.12  gives  a  vector  mcesure  of  the  inverse  detectability, 
a  measure  that  would  be  unwieldy  to  plot  or  even  record  for  very  many 
cases.  Application  to  a  specific  mission  or  circumstances  would 
usually  result  in  interest  in  only  one  or  two  variables  and  those 
would  be  selected  for  further  scrutiny.  In  the  absence  of  a  specific 
application,  Figure  7.2  plots  the  norm  of  the  position  error  produced 
divided  by  the  sos  increase. 

It  is  seen  that  there  are  wide  variations  in  the  detectability. 
For  example,  after  30  days  ar  would  produce  a  position  estimate  per¬ 
turbation  of  36OO  km  before  the  sos  rose  by  10C#.  For  the  same  sos 
change,  ax  would  produce  only  a  180  km  error.  Again  the  relative  re¬ 
lationship  between  r  and  \  effects  changes  as  more  data  is  added. 

c 


138 


E(SOS) 


i 


10s 


10" 


■THREE  ORTHOGONAL  ACCELERATIONS 

1  a  *  10  ^ 

2  *  10  kmX«c^ 

3  a  -  10  11  km^#c^ 

Y 

TWO  STATION  LOCATION  ERRORS 

4  r  *  10  m 

t  • 

I  5  X  ■»  10  m 
:  DOPPLER  ONLY 


[DATA  WEIGHT  AT  1  mm/i«c 
[  1  MINUTE  SAMPLE 


102 


0  5  10  15  20  25  30  35 


Figure  7.2.  Injection  Position  Error  Per  Unit  so*  Increase 
versus  Tracking  Time. 


Concluding  Remarks 

In  summary,  we  feel  that  the  considerations  discussed  in  this 
chapter  provide  at  least  the  genesis  of  a  more  formal  estimation 
validation  procedure  which  might  be  performed  as  follows: 

1.  Select  a  data  set  .aere  relevant  functions  of  the  d's  are 
adequately  low  for  all  the  consider  variables. 

2.  Process  the  data:  if  the  resulting  sos  conforms  well  to  that 
predicted  statistically  from  random  data  noise  considerations 
alone,  high  reliance  can  be  placed  on  the  orbit  obtained. 

3.  If  not,  inspect  the  data  rotated  by  the  appropriate  House¬ 
holder  which  also  separates  the  residuals  into  three  classes: 

I.  z  arising  from  state  variables,  consider 
x  variables,  and  data  noise 

II.  z  -  arising  from  consider  variables  and  data 
y  noise 

III.  z  -  arising  from  data  noise. 

These  classes  appear  in  the  rotated  data  vector  as  follows: 


Observe  the  two  residual  classes  z  and  z  .  If  the  cause  of  the  high 

y  z 

aos,  or  structured  residual  pattern,  is  largely  in  the  z  class,  this 

z 

means  that  what  is  causing  the  difficulty  has  not  been  modeled  even  in 
the  consider  parameters,  let  alone  the  solution  parameters.  In  this 
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case  there  is  little  that  can  be  done  other  than  reaching  the  unsatis- 

Tying  conclusion  that  the  difficulty  is  at  least  in  part  arising  from 

considerations  not  even  thought  of. 

4.  If  the  high  residuals  are  largely  contained  in  the  z  class, 

¥ 

the  consider  parameters  span  the  difficulty.  Execute  the  test  algor¬ 
ithm  (Equation  7.16)  to  determine  if  one  or  perhaps  a  few  consider 
parameters  are  present  at  a  high  value.  Once  the  culprit (s)  has  been 
isolated,  add  it  to  the  estimation  list  and  re-perform  the  procedure. 


APPENDIX 


I.  Orbit  Determination 

Orbit  determination  is  the  process  of  using  imperfect  observa¬ 
tions  which  are  functionally  dependent  on  the  state  of  the  observed 

body  to  form  an  estimate  of  the  body's  trajectory  which  is  best  in 

.  ’  -  •  "  i  ’  *'  k  •  ,  ! 

some  defined  sense.  We  shall  be  speaking  of  trajectories  subject 

primarily  to  gravitational  force  fields  and  in  particular,  space¬ 
craft  trajectories .  Further,  we  shall  be  concerned  primarily  with 

'  V'.’ 

earth-based  observations  of  the  spacecraft;  the  most  useful  of  these 
are  tracking-station  to  spacecraft  range  and  range  rate.  They  will 

be  discussed  in  more  detail  in  the  following  sub-section. 

.jvv.  ..•••  -/'•  ' 

'  V  4  ..  ■  '  ■ 

The  basic  process  of  using  the  measurements  to  form  the  esti- 

■  *_'< ..  *  ,  ’  *  '  ,  ,  ’v  ,  -  "'v. 

mate  is  shown  schematically  in  Figure  A.l.  Here,  our  point  of  view 
is  to  obtain  an  estimate  of  the  initial  conditions,  or  state,  and 
then  through  integration  of  these  estimates  form  the  complete  es- 

•'  *  ..  '  *  * .  jf  '  '  '  *  ■{ 

timated  trajectory.  It  is  often  convenient,  though  by  no  means  a 
requirement  j  to  visualize  the  initial  conditions  as  the  position  and 

,  -r.  .  ’  n 

, v  ••  ...  •  .•  *  I*-.  ’  -  .  J 

velocity  of  the  spacecraft  at  that  epoch.  More  generally,  the 

initial  conditions  will  Include  both  the  spacecraft  state  and  other 

parameters  that  affect  either  the  propogation  of  spacecraft  state 

and/or  the  measurements  themselves.  To  be  more  explicit  we  define: 

*  A  total  parameter  a  T  x(0) 
vector  being  g 
estimated 

8 

.  Preceding  page  blank 

»3 


INITIAL 

GUESS^  START  AGAIN 


Figure  A.l,  Functional  Diagram  of  Orbit  Determination  Process. 
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where  x(0)  is  the  spacecraft  state  (position  and  velocity)  a  six 
vector 

v  are  parameters  which  influence  the  propogation  of  spacecraft 
state  (e.g.  mass  of  the  earth,  solar  pressure) 
s  are  parameters  which  affect  only  the  data  (e.g.  station  lo¬ 
cations  ,  data  biases) 

It  is  important  to  restrict  the  members  of  q  to  only  those  para¬ 
meters  where  the  uncertainty  in  their  true  values  significantly 
affect  the  problem  at  hand.  Parameters  that  are  relatively  well 
known  or  whose  involvement  in  the  problem  is  minimal  are  treated  as 

’I  ’  *  ,  ■/ 

perfectly  known  and  do  not  appear  in  q  (e.g.  mass  of  Pluto  for  a 
Mars  Orbiter  problem). 

To  start  the  process,  an  initial  guess  of  q  must  be  furnished. 
This  is  used  to  integrate  the  non-linear  equations  of  motion  and  pro¬ 
duce  a  nominal  probe  ephemeris.  Then,  a  linearization  about  this 
nominal  takes  place  to  calculate  the  effects  of  small  variations  in 
q.  That  is,  the  system  of  equations 

t(t,to)  ■  A$(t,t0) 

where 


A 


2* 

A* 

ax 

av 

ii 

&X 

av 

is  evaluated  about  the  nominal 


in  a  straightforward  manner.  Since  the  v's  are  formulated  as  con¬ 
stant  parameters  of  the  system,  the  lower  partitions  of  A  are  zero 
and  the  $  obtained  has  the  special  form  . 
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The  probe  trajectory  is  also  used  to  calculate  "estimated" 
observables.  This  is  done  by  simply  asKing  the  physically  meaning¬ 
ful  question,  "If  the  actual  trajectory  were  exactly  as  I  have 
guessed  it  to  be,  and  a  perfect  observation  were  taken,  what  would 
be  the  value  of  that  observation?"  The  estimated  observation  is 
differenced  with  that  actually  received,  producing  residuals  that 
are  fed  directly  into  the  estimation  process. 

In  addition  to  the  residuals,  we  have  need  of  the  partial  de¬ 
rivatives  of  the  data  with  respect  to  the  total  parameter  vector,  q. 
Both  the  U  and  V  matrices  already  defined  are  needed  in  this  final 
computation.  Let  the  data  received  at  time  t  be 
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The  v  enters  both  explicitly  in  h  and  as  an  argument  of  x(t) 
since  certain  parameters  can  affect  both  the  propogation  of  the 
spacecraft  state  and  the  data  directly.  These  are  unusual  however 
and  do  not  arise  often;  perhaps  the  best  example  is  the  ephemeris  of 
the  Earth.  Its  position  obviously  affects  the  trajectory,  but  in 
addition,  a  change  in  the  Earth's  position  moves  the  tracking  sta¬ 
tion  and  contributes  directly  to  a  change  in  the  data. 

Together,  the  partiala  and  the  residuals  are  used  in  the  esti¬ 
mation  process  to  form  the  estimate  of  q,  "q.  Along  with  this  esti¬ 
mate  will  also  be  statistics  on  the  errors  in  the  estimate  in  the 
form  of  a  covariance  matrix,  for  the  second  moments  (the  first 
moments,  generally,  are  zero).  Usually,  the  quantity  of  prime  in¬ 
terest  is  not  the  initial  conditions,  but  the  spacecraft  trajectory 
near  or  at  a  distant  target  body.  If  the  later  is  the  case,  the 
estimated  Q  is  used  to  produce  a  new  complete  trajectory.  In  addi¬ 
tion,  the  statistics  on  the  state  can  also  be  propogated  by 


[u !  v]  a  fui  vlT 
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where  U  and  V  are  understood  to  have  arguments  t,  tQ. 
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II.  Earth- Based  Radio  Tracking  Data 
Definitions  and  Characteristics 

We  include  here  a  brief  description  of  the  physical  system  em¬ 
ployed  by  the  Deep  Space  Network  (DSN)  of  the  Jet  Propulsion  Labora¬ 
tory  and  a  measure  of  the  inherent  accuracy  of  that  system. 

Figure  A. 2  shows  the  basic  building  blocks  for  a  coherent  range 
and  doppler  system.  This  figure  could  have  been  drawn  as  separate 
diagrams  and,  indeed,  separate  systems  could  have  been  mechanized  for 
each  of  the  data  types.  It  is  useful,  however,  to  emphasize  not  only 
the  common  elements  of  the  hardware  system,  but  the  conceptual  rela¬ 
tionship  between  the  two  data  types.  Basically,  the  range  system 
measures  the  roundtrip  light  time  between  the  probe  and  station, 
while  the  doppler  system  measures  the  time  rate  of  change  of  this 
quantity.  The  latter  is  accomplished  by  sending  a  signal  having  a 
stable  frequency  to  the  spacecraft.  The  signal  received  at  the 
spacecraft  is  re-broadcast  to  the  station  and  is  received  with  a 
frequency  shift  proportional  to  the  relative  range-rate.  The. 
doppler  tone  produced  is  simply  integrated  by  means  of  an  electronic 
counter  which  results  in  the  Integral  of  the  range-rate. 

At  the  same  time,  the  carrier  is  phase  modulated  with  the  rang¬ 
ing  code.  The  code  broadcast  can  have  a  variety  of  mathematical 
structures  depending  on  the  formulation  used.  A  particularly  in¬ 
teresting  code,  however,  one  currently  being  used  by  the  DSN,  is  the 
so-called  multi-component  "pseudorandom  code."  Although  this  par- 
ticluar  code  is  generated  det''rr’irintica?ly  and  has  a  finite  length 

Preceding  page  blank 


Figure  A. 2.  Coherent  Range,  Doppler  Ground  Tracking  System. 
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before  repetition,  it  takes  on  many  of  the  characteristics  of  a  purely 
random  sequence  of  l's  and  O'a.  When  this  signal  is  ~e  turned  to  the 
station,  it  has  shifted  in  time  corresponding  to  the  instantaneous 
round-trip  light  time  involved.  The  signal  is  then  auto-correlated 
with  the  code  being  generated.  The  generated  code  is  also  shifted 
in  time  as  shown  in  Figure  A.  2.  Here,  the  properties  of  the  pseudo¬ 
random  sequence  are  used,  and  the  correlation  computed  will  be  uni¬ 
formly  low,  unless  the  generated  code  is  shifted  in  time  by  the  same 
amount  that  the  received  code  is  shifted  by  its  round-trip  transit 
time.  By  carefully  constructir^  the  code  from  subcodes  of  gradually 
increasing  lengths,  this  search  for  maximum  correlation  can  be  done 
efficiently  and  need  not  require  a  trial  correlation  with  every  poss¬ 
ible  shift  time.  The  amount  of  shift  required  to  produce  this  maxi¬ 
mum  correlation  becomes  the  raw  data  type  proportional  to  the  space¬ 
craft  range.  Therefore,  the  counted  doppler  is  seen  as  the  Integral 
of  the  range-rate  while  the  ranging  system  supplies  the  arbitrary 
constant  of  integration.  These  are  the  raw  data  types  actually  used. 

This  system,  particularly  the  doppler  portion,  is  relatively 
unsusceptible  to  instrument  error.  The  resulting  high  inherent  data 
accuracy  has  deferred  the  need  for  on-board  tracking  systems,  and  has 
propelled  the  Earth-based  system  as  the  prime  means  of  orbit  determin¬ 
ation  far  beyond  the  limits  originally  envisioned  when  planetary  ex¬ 
ploration  was  first  undertaken. 

In  order  to  establish  a  means  of  comparison,  consider  tracking  a 
spacecraft  having  a  typical  Earth-relative  velocity  of  10  km/sec. 

Since  the  raw  data  type  is  ran^e  change,  and  this  is  measured 
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continuously  over  periods  up  to  12  hr  (s  single  pass),  the  measure¬ 
ment  will  be  on  the  order  of  $00,000  km. 

Perhaps  the  most  visible  noises  impressed  on  the  data  are  re¬ 
ceiver  jitter  and  rather  high-frequency  variations  arising  from  the 
master  oscillator  instability.  These  errors  are  visible  because  of 
the  high-frequency  content,  and  will  be  in  evidence  after  even  a 
rudimentary  smoothing  of  the  data.  Their  excursions  are  small,  how¬ 
ever,  not  amounting  to  more  than  a  few  centimeters,  or  an  upper 
bound  of  one  part  in  10^  of  the  basic  measurement.  In  addition,  be¬ 
cause  the  noise  is  high  frequency  tends  to  make  it  easily  filtered 
out  in  any  orbit  determination  process,  and,  therefore,  it  does  not 
seriously  affect  the  estimates. 

Low-frequency  errors,  from  a  range-rate  bias  to  frequencies 
approximately  10  times  the  frequency  of  the  Earth's  rotation,  are  of 

:  •  !  •  •  •  •  r  v  •  ,  , 

greater  Importance.  Fortunately,  the  Instrument  is  not  particularly 
susceptible  to  these  errors.  Consider  again  Figure  A. 2  and  note  that 
a  bias  in  the  oscillator  frequency  will  cause  a  shift  both  in  the 
transmitted  and  received  frequencies,  quantities  that  are  differenced 
to  produce  the  doppler  tone.  The  effect  of  the  bias  is  nearly  can¬ 
celled  by  this  operation.  If,  on  the  other  hand,  a  change  in  the 
frequency  occurs,  the  doppler  will  be  affected  directly  by  the  magni¬ 
tude  of  that  change.  The  shift  in  the  doppler  tone  is  integrated  to 
produce  a  range-change  error.  A  cancellation  occurs  later,  however, 
when  the  frequency  change  reappears  in  the  received  signal;  this  time 
it  enters  the  data  in  the  opposite  sense.  Therefore,  an  oscillator 
bias  is  of  no  consequence;  it  is  the  variations  ir.  the  transmitted 
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frequency  that  produce  first-order  effects.  Current  rubidium  s tan- 

12 

dard  performance  of  drifts  up  to  5  parts  in  10  over  a  pass  will 
produce  errors  of  80  cm  for  a  20-minute  round-trip  light  time  (as 
shown  in  Chapter  II). 

The  approximate  doppler  equations  shown  in  Figure  A.  2  give  the 
impression  that  the  uncertainty  in  the  speed  of  light,  c,  would  give 
rise  to  a  bie  •  in  the  range-rate  measurement.  This  arises,  however, 
only  because  c  enters  from  a  desire  to  express  velocities  in  metric 
units;  it  is  not  a  requirement.  All  measurements  and  ephemerides  are 
fundamentally  taken  and  calculated  in  light-seconds,  the  distance 
light  travels  in  unit  time.  Therefore,  c  is  a  defined  constant  for 
convenience  only;  it  could  as  well  be  unity. 

There  are  a  variety  of  sources  giving  rise  to  errors  in  the  re¬ 
mainder  of  the  low-frequency  region.  In  total,  they  add  to  a  stan- 
dard  range-change  error  of  approximately  1.5  m  during  a  single  dop¬ 
pler  tracking  pass,  or  using  the  original  yardstick  of  500,000  km, 

3  parts  in  109. 

The  current  planetary  ranging  system  is  capable  of  measuring  the 
station-spacecraft  distance  to  within  a  standard  error  of  approximat¬ 
ely  15  m  bias  plus  a  5  m  random  component.  These  numbers  justify  the 
earlier  statement  that  the  ranging  data  supply  the  constant  of  range- 
rate  integration;  i.e.,  on  the  basis  of  a  single  pass,  the  doppler 
gives  a  considerably  more  accurate  measurement  of  range  change  than 
does  the  difference  of  two  range  measurements.  Taking  a  longer  view, 
the  situation  is  reversed  since  doppler  errors  villcontinue  to  accu¬ 
mulate  but  the  difference  betv*»er.  two  rnr’*o  measurements  ever,  after 
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several  days  will  still  be  accurate,  to,  say,  'JZTx  5  m.  To  capitalize 
on  this  feature,  one  oust  pay  careful  attention  to  maintaining  the 
calibration  of  the  ranging  system  to  ensure  that  the  bias  is  liter¬ 
ally  that,  and  not  a  slowly  varying  quantity. 


III.  The  Estimation  Process 


We  now  turn  our  attention  to  the  estimation  process  itself. 

Let  the  data  be  represented  as  a  linear  function  of  the  state  (actu¬ 
ally,  of  course,  the  perturbation  from  nominal  in  the  state)  plus 
additive  noise,  i.e. 

z  ■  Hx  +  c  (A.l) 


where  z  is  an  m  vector 

x  is  the  (unknown)  state,  an  n  vector* 
and  H  is  an  mxn  matrix  of  partial  derivatives 
<  is  an  m  vector  of  the  measurement  error. 

It  is  assumed  that  ■  0  and  e£*  cTJ  ■  Af  where  E  is  the  expected 
value  operator;  superscript  T  denotes  the  transpose.  Our  objective 
la  to  find  a  linear,  un-biased  minimum  variance  estimate  of  x,  x. 


Or  stated  precisely 

i.  x  ■  Bz  (linearity) 

il.  E[A]  -  x  (unbiased) 

iii.  a  -  Ef<x-*)(x-ft)Tl  is  a  minimum 
*  L 

Substituting  (4.1)  into  i  and  employing  ii 

x  -  B(Hx  +  s)  (A.2) 


Epcj  »  BHx  ■  x  if  and  only  if 

BH  -  I  (A.  3) 


— 

*It  is  convenient  to  now  drop  the  distinction  between  x,  v,  and  s  and 
treat  all  parameters  in  a  unified  fashion.  Thus,  x  here  is  the  same 
as  gq  in  Section  I  end  is  used  to  conform  with  more  universal 
notation. 


155 


Substituting  (A.P)  and  (A. 3)  into  iii  gives 

A*  m  e[b«  «V]  -  Ba#BT  (A. 4) 

Equation  (A.4)  must  be  nininized  with  respect  to  B  but  observing  the 
constraint,  (A. 3).  Thus,  adjoining  this  constraint  with  &  matrix 
Lagrange  multiplier  to  (A.4) 

A*  -  BAcBT  +  (BH-l)x  +  XW  -  I)  (A. 5) 

In  order  for  (4)  to  be  a  minimum,  the  variation  of  (5)  must  be  zero. 

6AX  “  *b<a«bT  +  Hx)  +  +  xVjiB1  -  0  (A. 6) 

Evidently 

BAf  ♦  xV  ■  0 


or 


and  since  BH  -  I 


thus 

x  -  (A.7) 

and 

Ax  -  (A. 8) 

These  are  the  familiar  weighted  least-squares  or  batch  minimum  vari¬ 
ance  estimation  formulae.  They  hove  been  in  use  since  Gauss  with 
very  little  modification  and  arc  still  the  mainstay  of  many  sophis- 

'  i 

ticated  orbit  determination  problems  in  use  throughout  the  country'. 

With  very  little  additional  work  we  can  introduce  the  concept  of 
apriori  information  on  the  state  prior  to  data  taking. 


Suppose  we  have  an  initial  guess  as  to  the  value  of  x,  in  the 


form 


x  “  x+«x 

where  e„  is  mean  zero  with  covariance  a  . 
x  "ap 

This  can  be  considered  as  simply  additional  data  related  to  x  through 
an  identity  transformation.  Since  usually  (though  not  necessarily) 
the  initial  guess  at  the  trajectory  will  correspond  to  the  apriori 
data,  in  our  notation  that  data  becomes  zero}  i.e. 


(A.9) 


”  0  " 

*  I  " 

e 

X 

B 

x  + 

z 

H 

« 

m  mm 

mm  — 

Substituting  (A.9)  into  (A.7)  and  (A. 8)  yields 

-1 


m 


rsi 

'ap 


0 

m  mmm 
-1 


I 

■  M  m 

H 


.  r*Ip  0  ‘ 

r  I !  HT1  I - »— 

L  '  0  a;1 


and 

tT.-l,,\-l 


Ax  *  <Alp  ♦  »V>) 


(A.7) 1 
(A.8)' 


Sequential  or  Kalnan-Bucy  Filtering 

A  version  of  the  Gauss  Markov  Theorem  (Reference 2k)  states  that 

if  two  mean  zero  random  variables,  x  and  z,  have  covariances  a  and 

"ap 

•  '  ’  .  •  - 

A  respectively  and  cross  covariance  a  ,  then  the  minimum  variance, 

Z  aZ 

'  •  •  ?  .  *  —  * 

linear,  unbiased  estimate  of  x  given  z  is: 

A  -1 

x  -  - - 


»  fm  "N. 


^A»10) 
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(A. 11) 


and  the  resulting  covariance  in  the  error  £  is  given  by 

-IT 

Ax  *  A*p  *  ^xz^zz  AxZ 
If  z  is  again  considered  the  data  in  accordance  with  (1),  and  x  has 
prior  uncertainty  A  t  then 

^ G  *  ,  O  4  •  *  f  jh ri  »  r  •  •  ■  ■  *  \ 

Axz  "  AapH 

A*.  ■  («  A,/  ♦  A.) 

Thus (A.  10)  and(A.ll)becooe: 

*  “  AapHT(«  AapH^  ♦  Ac)”lz  (A.  10) ' 


Ax  ’  A.p  -  V'T(H  A.PhT  +  A,)*1*  App 


(A.  11)' 


This  in  essence ,  is  the  Kalman-Bucy  filter  in  discrete  form  for  the 
very  simplest  of  circumstances  (Reference 25).  We  do  not  present  in 
detail  here  the  equations  for  the  full  recursive  algorithm  including 

L  J  .  i  ‘  $ 

mapping  between  data  times  because  they  are  lengthy  and  only  serve  to 
obscure  the  elegance  of  the  principal  results  -  all  of  the  basic 

ideas  are  contained  in  these  two  equations.  In  words,  the  procedure 

. 

becomes: 

If  (A. 10)’  and  (A.ll)'  represent  the  results  from  the  first 
data  point  (or  points)  and  it  is  desired  to  process  more  data  - 

l)  replace  with  corrupt  ^  by  the  addition  of  a  positive 

definite,  symmetric  matrix  if  random  process  noise  is  enter¬ 


ing  the  system 

2)  subtract  from  the  new  data  H'ft,  which  is  the  expected  value 
of  the  new  data  (rf^is  TTflte  new^arafrol  derAuqt^^matrix)  ^ 

3)  re-execute  (10)’,  (11)';  the  resultin';  r  is  the  ".mount  the 
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estimate  should  be  changed  from  the  previous  value. 

Both  the  least-squares  and  sequential  methods  are  unbiased,  mini¬ 
mum  variance  techniques.  If  this  is  so,  it  must  be  possible  to  show 
the  equivalence  between  (A.7)’ ,  (A. 8)',  and  (A.10)',  (A. 11)'.  We  be¬ 
gin  with  an  apparent  diversion.  Suppose  we  have  two  matrices,  par¬ 
titioned  commensurably,  one  the  inverse  of  the  other 


*  »  1  r  «  8*1  fx  ;  o 

c  ®  J  L Y  8 J  "  L‘°Tr 

Then 


Aa  +  By  -  I  (i) 

Ca  +  Dy  «  0  (ii) 

First,  solving  (ii)  for  y,  and  using  in  (i),  we  obtain 

Y  -  -D*’1Ca 

a  -  (A  -  BD^C)”1  (iii) 

Y  -  D”*C(A  -  BD’^C)’1  (iv) 


Second,  solving  (i)  for  a  and  using  in  (ii) 
a  -  A"1  -  A_1By 

a  -  A’1  -  A'^CA^B  -  D)"1CA*1  (v) 

Y  -  (CA~1B-D)’,1CA”1  (vi) 

Equating  (iii),  (v)  and(Lv),  (vi)  ve  have  finally 

(A  -  BD^C)”1  -  A-1  -  A"1B(CA"1B-n)_1CA‘1 
(CA"1B-D)_1CA‘1  -  -D-1C(A  -  BD^C)"1 
Now  for  (A.  12)  if  we  let 


D 


(A.12) 
(A. 13) 


B 


H 


A  -  Ae 

Then 

(HTa;1h  *  A^'^V  -  A./U,  +  HA,/)-11 

showing  the  equivalence  of  the  two  filter  equations. 
Finally  if  for  (A. 13)  we  let 


C  -  H 

Then 

<A;J  ♦  X^1")'1  -  Aap  -  Aa/(HAHT  ♦  At)'Xp 
showing  the  equivalence  of  the  two  covariance  equations. 


Scuare-Root  Filtering 

The  basic  motivation  of  snuare-root  filtering,  which  can  be 
used  either  with  the  least-squares  or  sequential  techniques  is  that 
by  operating  with  either  the  square  root  of  the  covariance  or  infor¬ 
mation  (inverse  covariance)  matrices,  numerical  integrity  is  improved. 
This  is  best  illustrated  by  a  specific  example.  Suppose,  we  wished 
to  solve  the  following  equation  for  x 

z  ■  Hx 


on  a  digital  computer  in  single  precision  (8  decimal  places).  To  do 

this  in  the  usual  way  would  involve  taking  the  determinant, 

det  H  ■  1  +  c  -  1 
-6 

If  c  were,  say  «  10  ,  the  calculated  determinant  would  have  two 
place  accuracy  and  a  reasonable  answer  would  result.’  If,  however,  we 
were  to  solve  this  problem  in  a  least-square  sense  which  theoretically 
yields  the  sane  answer,  i.e.,  if  we  were  to  compute 
x  -  (HTH)-1  HT  z, 

T 

we  would  have  to  essentially  find  the  determinant  of  H  H,  which  is 
readily  seen  to  be 

det  HTH  «  U  +  he  +  2c2  -  4  -  Uc  -  t 

If  e  were  10-^,  this  would  compute  as  zero,  and  the  matrix 
could  not  be  inverted.  Since  for  accurate  orbit  determination  it  is 
necessary  to  include  as  estimated  variables  all  parameters  that  con¬ 
tribute  information  to  the  data,  the  problem  of  a  nearly  rank  defic- 

T 

ient  H  is  chronic.  Thus  the  motivation  is  high  to  avoid  the  H  H  op¬ 
eration  and  to  find  ways  of  dealing  with  H  directly  even  (particularly) 
when  it  is  non-square.  We  present  the  outline  of  one  such  technique 
here  that  employs  a  generalization  of  the  Householder  trans format ions 
originally  used  for  inverting  square  matrices  by  first  triangulariz- 
ing  and  then  inverting  by  backwards  substitution.* 

We  begin  with  the  least-squares  data  model 

z  -  Hx  +  t  (A.  1*») 

•Definitions:  Dy  th?  rquare-root  of  a  square  matrix  Q,  we  simply 
moan  th"  (non-uni cu*1)  Tatiix  such  that  Q2  qzt  ■  Q.  An  upper  (lover 
triangular  matrix  has  all  zero  elements  bclow(above)  the  main 
diagonal. 
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where  the  sane  definitions  apply  as  before. 

If  (A.l4)  is  pre-multiplied  by  A^  and 

A  -i 
a  -  A*  c 

a"^  Z  -  X  +  n 


(A.15) 


where  now 


E[n]  =  0,  E  £n  nTJ 


V/e  seek  an  estimate  of  the  state  x  and  the  noise,  n,  such  that  ft  is 


minimum  subject  to  Equation  (A.15),  i.e. 


+  n  -  A‘*  * 


(A.16) 


It  is  possible,  through  an  orthonormal  transformation  P  to  find  a  new 
basis  such  that  the  column  vectors  of  A~^  H  are  contained  within  the 
sub-space  formed  by  the  first  n  vectors  of  the  new  basis,  i.e. 


f  ^  H  - 


m-n  0 


Equation  (A.l6)  becomes 


—u  ♦  — 


(A. 17)* 


8ince  P  is  orthonormal,  «  llft|l  and  the  problem  is  still  to 

ft^  r*  -jj. 

satisfy  (A.17)  and  minimize  the  J.  This  solution  is  simply 

IK} 


.  Thin  solution  is  simply 


*Tho  conditions  on  P  thus  far  imposed  mnke  ?  non-unique.  We  can 
facilitate  the  later  inversion  of  H  and  cake  P  unique  if,  in  addi¬ 
tion  we  require  that  R  be  upper  triangular 


K  -  0 


4  ■  R'll! 


To  show  equivalence  with  the  least  square  formulation  just  discussed  we 


observe: 


1.  The  solution  is  linear 


2.  it  is  unbiasedy  since 


k  ■  R"1*^ 


R-1Rx  +  R_1  Pn 


eT*J  -  lx  +  R-1P  sTnl  -  x 

3*  the  covariance  of  the  estimate  is  still  the  sane. 
x-k  -  R_1Pn 

£f(x-ft)(x-i)Tl  A  .  R_1P  Efn  n*J  fV1*-  R'V1* 

T 

since  P  P  -  I  froa  the  consequence  of  the  orthononeeUty. 


But 

[*TI  0  ]  [7]  •  *  p  * 


R'V1  -  (H1  a"1")"1 


which  is  the  sane  as  the  weighted  least  square  covari- 
anct  and,  by  the  proof  of  the  previous  sub-section,  as 
the  sequential  foneula. 
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